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ABSTRACT 


As continued efforts are made to reduce the radar cross sections of aircraft and 
ships, designs are first modeled with computers and then tested in the lab. In the far-field 
of these tested objects, actual sources of high reflectivity or “Hot Spots” on the tested 
objects can be isolated to within only one half the wavelength of the electromagnetic 
wave used for testing. Ideally, a probe could measure fields on the surface of the object 
being tested to completely isolate the source of the hot spot. Unfortunately, the presence 
of the probe on the surface of the object will disturb the very fields it is attempting to 
measure. Probe measurements made in the near field, close to but not on the object, can 
be designed to reduce the influence of the probe while providing accurate field data. The 
data thus measured, while not able to determine the source location perfectly, can be 
used to localize a source to less than one half wavelength, the far-field diffraction limit. 
This thesis tests a technique for back propagating computer generated near field 
measurements of an axisymmetric field source to determine the fields closer to the 
source. Several cases are examined that test the accuracy and resolving capability of the 
technique. 
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L INTRODUCTION 


The harder it is to detect a platform, the easier it is for a platform to escape 
observation and targeting. This is the driving force behind research and engineering in 
the area of “stealth.” Stealth refers not only to the oddly shaped airframe of the F-117, 
but to a wide realm of techniques designed to reduce the many facets of a platform’s 
“signature.” The signature of a platform consists of acoustic, thermal, and 
electromagnetic sources that are exploited in the detection of that platform. As stealth 
technology has become more publicized, various efforts for coimtering it have also been 
under exploration. Advances made in overcoming stealth directly imply that continued 
efforts are required to further reduce signatures. 

A. GOAL OF RESEARCH 

In the areas of Radar Cross Section (RCS) and electromagnetic sources, the 
general trend has been to attempt back-propagation of field measurements from the far 
field region of a platform down to the surface of the platform itself From this estimation 
of the fields on the platform, the distribution, localization, and characterization of source 
cuirents can be approximated. Unfortunately, sources can only be distinguished 
separately if they are greater than one half wavelength apart when measurements are 
taken from the far field. Figure 1 shows the relationship of three regions of interest to the 
surface of a source and the regions’ resolution limitations. 

Alternatively, since exact source localization would be ideal, field measxuements 
made directly on the platform by a perfect probe might yield the required information. 
Unfortunately, the presence of an actual probe will greatly disturb the fields and source 
distribution on the platform. 

This leaves the near field region of the platform, between the platform’s surface 
and the far field, in which to make measurements and to image the sources. For a probe 
in the near field, the measurable fields no longer constrain the source resolution to 
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structures greater than one half wavelength (as was the case for the far field) and allow 
for greater refinement. Additionally, the probe can be placed at a distance sufficiently far 
away from the surface to make the ejBfects of its presence managable. 


Multiple 
sources with 
separation d 



Figure 1. Far Field/Near Field/Surface Localization Limitations 

Exploring a technique for the back propagation of fields measured in the near 
field region at one distance from the source onto fields closer to the source is the primaiy 
goal of this thesis. 

B. BACKGROUND OF RESEARCH 

The analysis of near field phenomena that spurred the current research originated 
in the field of acoustics. For plate and cylinder like structures, a newly defined quantity 
called the supersonic acoustic intensity vector provided a tool for the location and 
characterization of acoustic sources from near field acoustic measurements. [Ref 1] 
From acoustics, the supersonic analysis of sound wave components was translated into 
the realm of electromagnetics. For cylindrical geometiy, the now superluminal modes 
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which satisfy |A:j| < ^ and have faster-than-li^t propagation along the cylindrical axis 

provide time-average power to the far field. Measurements taken in the near field were 
back-propagated using an optimal deconvolution filter to account for the subluminal 
modes. [Ref 2] This technique of using superluminal modes was then extended from 
separable geometries to non-separable geometries. Using the Finite Element Method to 
relate the boundary conditions of a measured field to the surface currents induced on an 
axisymmetric body by scattering, fields were back-propagated to estimate the source 
distribution. Resonance issues in the spherical test case complicated the overall 
effectiveness of that particular approach. [Ref 3] 

This thesis continues to explore methods of near field back-propagation that apply 
to general non-separable geometries from axisymmetric bodies. These are bodies that do 
not have explicit separable modes such as cylindrical or spherical wave functions. 
Unlike the previous work in non-separable geometries, this thesis attempts to reconstruct 
fields closer to the source and not strictly the source itself. Additionally, this thesis looks 
at source current distributions that are not induced by scattering. 

This study is organized as follows: Various cases of generated electromagnetic 
radiation source fields are discussed in Chapter H. Various techniques in calculating 
forward propagation fields and testing their convergence are examined in Chapter in. 
Various techniques in calculating back propagation of fields are discussed in Chapter IV. 
A summary and concluding discussion are presented in Chapter V. An appendix has 
been included with the fimdamental set of computer codes to allow for an indepth review 
and continuation of this line of research. 
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n. GENERATION OF ELECTROMAGNETIC SOURCE FIELDS AND THE 

FIELD TRANSFORMATION OPERATOR 

A. GENERATION OF THE EXACT FIELDS 

Before one can perform numerical experiments on the accuracy of back- 
propagated fields, one must generate the fields to be propagated and transform matrix 
that does the propagation operations. First we consider a single dipole located at the 
origin. Following the derivation [Ref. 4], a closed form solution to the radiation patterns 
any distance away from the source exists. These closed form solutions form the exact 
fields at various p values which are in turn used for two pxuposes. First, the exact fields 
are used as checks of the accuracy for various numerical techniques, and second, as the 
basis for determining the equivalent electric and magnetic currents required for 
numerical propagation. The dipole is placed at the origin with its length along the z-axis 
as shown in Figure 2. 
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The axisymmetric fields it produces will be the same at any point with the same p 
and z around the axis. Any combination of dipoles oriented along the z-axis will also 
produce axisymmetric fields. The closed form solutions are based on an analysis of a 
dipole with a sinusoidal current distribution. Figure 3 shows the dipole geometry used in 
the derivation. The current distribution on the dipole is sinusoidual such that 



Figure 3. Dipole Exact Field Derivation Geometry 

Ig = where z' e \¥h-h\. Figure 4 illustrates two cases. The first 

case shows the dipole with an overall length less then one half the wavelength associated 
with the current frequency. The second case shows the dipole with a length between one 
half and one wavelength. Phaser Fields may be determined through use of the vector 
potential, computed for line currents as [Ref 5] 

^jkR 

A(x,y,z) = -^\lg{x',y\z’)^-^dz' (1) 

with R = ^jx^ +y^ +(z-z')^ . (2) 
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Current Distrubutions for Several Frequencies with Dipole Length <= 7J2\ Dipole Length = 0.1 



Current Distrubutlons for Several Frequencies with Dipole Length > 7J2\ Dipole Length = 0.1 



Figure 4. Sinusodial Current Distribution 
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Shifting to cylindrical coordinates and substituting the sinusoidal current distribution of 
the dipole for TXx',y',z'), the integration may be readily performed. Using Maxwell’s 

equations, the magnetic field due to vector potential A{p,^,z) is given by 

H. = — (v X a). The resulting electric field is given by £ = —^— (v x h\ The fields 
/r/ ’ ja>s^ 


produced by a sinusoidal current distribution can then be derived as the following, 

U ^ +e-M, _(2cosi*y'*^l 

* j4^p ‘ 


~ DIPOLE 

~ DIPOLE 





-jkR, “ 

( 9 cos _ 

.1 ^ j 

1 P ) 

Rb 

vvFO A/# * J 


+—-- \2coskh)- 


R. 


R 






R 


(3) 

(4) 

(5) 


These three equations are fimctions of 

1) Dipole height, h, 

2nf 

2) Wavelength, k = -, 

c 

3) Maximum Current Magnitude, Idipole> the 

4) Observation point in p and z, but not (|) since the fields are axisymmetric. 
Using Hj, and E^, the fields tangential to a cylindrical surface of any radius p can be 
determined within machine roimd off accuracy. 


B. TANGENTIAL FIELD GENERATION FROM EXACT AXISYMMETRIC 
FIELDS ON A CYLINDRICAL SURFACE 


The fields on a cylindrical surface at p 2 can be determined in two ways. First, the 
fields can be calculated from the actual source axisymmetric current distribution. 
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Second, if the fields are exactly known on some closed surface that contains the actual 
sources and is inside the p 2 cylindrical surface, then the surface equivalence theorem can 
be used to transform the known fields into electric and magnetic currents on the closed 
surface. Since the source current distribution is axisymmetric, the inner surface is chosen 
to be an infinite cylindrical surface. Figure 5 shows the relative geometry of the two 
coaxial cylindrical surfaces used in determining the fields at P 2 . 



Figure 5. Coaxial Cylindrical Surfaces Geometry 

The exact fields from the dipole source have both E- and H-field components 
tangential to a cylindrical surface as shown in equations (3) and (5). Applying the 
surface equivalence theorem, the fields at pi, become the following equivalent surface 
currents. 
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(7) 


=Exh^= E^{p^,<^z^z X pj = E^^ = MA 


Using Jsurface and Msurface, the vector potentials A and F can be determined via the 
inhomogeneous vector potential equations. 



Substituting A and F into Maxwell’s equations yields the following component fields 
[Ref5]. 


Ha = 


Ea=\ 




NxA 






VxHa 


Ef=-—VxF 

s. 


Hf = 


—VxEf 




H = Ha+Hf 
E = Ea + Ef 


( 10 ) 

( 11 ) 

( 12 ) 

(13) 

(14) 


— - e-^ 

In order to determine and E, ultimately the Green’s kernel - will have to be 

R 

integrated. The following subsection sets up the problem for that integration. 
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1 . 


Setting Up the Problem 


The geometry for an arbitrary field point is shown in Figure 6. 



Figure 6. General Case for an Arbitrary Field Point (p 2 , ^> 2 , Z 2 ) 

The steps for setting up the problem for an arbitrary field point are shown below. 


p = cos^2^ + sin^Jji' 

^ = cos(90 + + sin(90 + (j>^y (16) 

p' = cosmic + sin (17) 

= cos(90 + + sin(90 + (18) 
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/ Ad A# 

r =/?!/? +v 
r = p^p + z^z 


R = r-r’ = p^p-p^p’ + (z 2 - Zi)f 

R = \p 2 -P\ cos(<4 - <^ 2 )]/^ + Loi(- sm(<4 - <f>^ + (zj - Zj )z 


J = /,z' = y,z 






(19) 

( 20 ) 
( 21 ) 
( 22 ) 

(23) 

(24) 

(25) 

(26) 


There does not appear to be a simple closed form solution to the Green’s kernel 
integration over the infinite cylindrical surface required to find A and F. However, 
there is a straightforward numerical technique for solving the integration. First one 
divides the cylinder at pi into many small bands or rings that have a constant current 
magnitude on the ring as shown in Figure 7. Then the ring is subdivided into many 
smaller squares such that both electric and magnetic equivalent currents are constant on 
the square [Figure 8]. 


12 





Az corresponds to the width 
of the thin ring such that 

Az ^ ^ Az 

z,-< z, < z, + — 

' 2 ' 2 


Figure 7. Physical Meaning of Az 



Figure 8. Subdivision of the Ring Showing Constant Equivalent Currents 
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The influence of the currents in each square on a ring to the field at pz will be 
summed together to give the total influence of the ring. Since the fields at all p are 
axially symmetric, the influence of whole rings can be used when setting up the 
mathematical (matrix) representation of the system, as opposed to the individual 
influences of each small square on each ring separately. One sums the respective 
contributions to an individual point on the pz cylindrical surface made by each ring along 
the length of the inner cylinder to get an estimate of the actual field at that point on pz- 
This is repeated for each point on P 2 in order to obtain the estimated fields on the entire 
Pz cylinder. 

Taking advantage of the axial symmetry of the fields, the field point to be 
estimated is placed on the x-z plane, where y=0 and ^ 2 "^ [Figure 9]. [Ref 6] When the 
field point is on the x-z plane, the problem set up is simplified. 


il 

(27) 

li 

(28) 

f = z 

(29) 

M = (m^ sin(- ))p + (m^ cos 

(30) 

^=lP2-Pi cos^ijo + !/>,(- sin^,)]j^ + [z2 - Zi]z 

(31) 

R = ^pI + pI - 2piP2 cos(l\ + (z2 - zj 

(32) 
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Figure 9. Field Point (xa, 0, Z 2 ) and the Cylinder Segment 


By using a ring of width Azi, A is calculated by the summation of AAz segments. 

-M, 


AJ = -^ 
" 4;r 




r=*-ao_„ 


(33) 


Therefore, the segment of Az==(AA 2 )i is as follows. 

00 « 

In the limit of Az -> dz, AAz -> dAz, and = A ■ 


(34) 


i=-oo 


" 4^ j " p ‘ 


R 


(35) 


Following a similar set up F and AF are found to be the following. 

-JkR « ® ^ L _ \ 

.-‘'■fez 


^ ^ Pid(l),Az (36) 
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( 37 ) 


— s,Pxt^ A A 

4;r 4r 

Once the vector potentials have been broken into segments for integration, the 
determination for the electric and magnetic field segments to be summed over the inner 
and outer cylinders is the next set of steps. 

2. The Equivalent Electric Current Case and the Vector A Potential 


Using A/4^as found in equation (35), can be determined by the following 
equation. 


1 


1 


( 


AH. =—V><A^ = —Vx 


P,Pi^ f 
4;r 


■JkR 




(38) 


Since the curl operation is with regard to the coordinate system of the observer (p 2 , (1)2, Z 2 ) 
and the integration is with respect to <|)i, then the curl can be brought inside the 
integration as shown below. 


AH, 


Pi^ J 


4;r 

Pi^ f 

4;r 


-JkR ^ 


Vx m = 




R 


7t [ 

Y 

p 

i\ 

V- x / f 


j 1 

-fT 

il ^ J J 

L 


-JkR 


(VxJ,z) 




(39) 


Since JJ is a function dependent only on pi and Zi and the curl operation is with regard 
to the observer’s coordinate system (p 2 , (j)©, Z 2 ), then V xj^z=^0. One will encounter 
many varieties of the operation shown below. [Ref. 6] 


g-JkR 


R 




,»+2 




+1 


e-M (-\ -jk 


This current case yields the following for n-1, V —— = -^+ 

K \K K 


(40) 


-JkRR^ 
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Therefore, one can say the following. 




4;r 

Pl^z 


R 






R^ R^ 


sin^e ^’^d(l>^ ] 


(45) 


Since R = pf- 2pxp2c^^<k + ^J , it is an even function in <l)i over 

the range -ti to +Ji. All functions of only R are even. All functions of R and cos(j>i are 
also even. But, since sin<l)i is an odd function, even functions of R and odd functions of 
sin(j)i will be odd over the interval and will integrate to zero. Therefore, (AHa)p = 0. 
Because one will encounter various versions of the functions of R and combinations of 
functions of R and cos(j)i that are to be integrated, one defines Un and Vq as well as Wn 
(which will be used in the next section), as follows. jRef. 6] 


^ -jkR 

(46) 

St 


(47) 


W„= \{ cos^i f 

(48) 


Un, Vn, and are solved numerically by dividing the ring into many segments [Figure 
5] and summing their individual components. The Matlab programs INTPHI*.M does 
this numerical integration step. Since and are uniform around the ring and the 
differences arise from the physical distance to each segment and cos<l)i, then the segments 
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can be simuned around the ring. The ring will be the basic element for determining 
fields as opposed to each A(j)ixAz segment. Overall, this reduces the computation time 
and simplifies the problem set up. 

Substituting Un, V„, and with applicable values for n yields the following. 


(Aff pfy,+jkvJi 


(49) 


3. The Equivalent Magnetic Current Case and Vector F Potential 


The mathematical details of the equivalent magnetic current case for Ha and the 
equivalent electric current for Ep are very similar in form. 


AEa =-—VxAF 


(50) 


In AA, there was only a z component in that contributed to AH^. In AF there are 
two components of : a ^component and a ^component. These two components add 


some complexity to the previous mathematical derivation of AH^. 

d(f>i 


At: 

-;r1 



r 

r 



f 

<1 

X 




J 

-7t 


K 

1 R J) 





xM, 


— _ p,Az }(-l^-JkY 


AE, 


At: 


J 


rJkR 




R 


dA 

^x'M)d<l>, 


(51) 


(52) 


(53) 


Looking at the Rx term gives the following result. 

-2iXcos^,)^-( 22 -2iXsin^i)i^ + (F2C0S^i (54) 

Substituting RxM^ into AEp and breaking it into components yields the following 
result. [Ref 6] 
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(55) 


(56) 


(57) 

(58) 

(59) 

4. Cross Field Influence Term 

In Section n.B.2 “The Equivalent Electric Current Case and the Vector A 
Potential” set fourth above, AH^ at P 2 was calculated from A4, which was calculated 
from J which, in turn, was determined from a knowledge of the H-field at pi. 
Therefore, AH^ is the part of the total Aff at p 2 determined strictly from the H-field at 
Pi- Similarly in Section in.B.3 “The Equivalent Magnetic Current Case and Vector 
F Potential,” AE^ at p 2 was calculated from AF , which was calculated from M, 

which was determined from knowledge of the E-field at pi. Therefore, AEp. is the part 
of the total AF at P 2 determined strictly from the E-field at pj. But these are only parts 
of their respective total fields. Since AH jotal = AH^+ AHp and AE^ool =AEp+AE^ , 
the influence on AH from the E-field at pi and the influence on AE from the H-field at pi 
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must be determined in order to accurately predict the fields at p 2 . Three ways to 
determine these field components were calculated. 

0. Method A: Maxwell’s Equations Applied Directly to Calculated 
AH^ and AEp Fields 

This method, the most strai^tforward, takes the three previously 
calculated equations of AH^ and AEp and computes AE^ and AH^ from them. Using 
Maxwell’s equations and substituting the calculated AH^, one can state the following. 


AE.=-^VxAH = ^—VxAHJ ( 60 ) 

jcos, jcos. 



Since all the representations contain partial derivatives of Un and Vn with respect to pz 
and Z 2 , it is useful to calculate them in terms of other Un, Vn, and Wn functions as shown 
below. [Ref 6] 

^C/„ = MnV„^2 + ( 64 ) 

dp^ 

dp^ 


( 65 ) 



( 66 ) 


dz2 

■^V„= -fe - 2i X«^„+2 + ) 

azj 

Substituting the partial derivatives of Un and Vn into AE^ 
following results. 

(AEj, = ^SjkU,-k%\ 

^ jcos.An: 

- p)^V,+ 2 jkV,-k%\} 


(AEj. = + jkU,\- .^[Kj + JkV^] 

jcos.Ak p2 

+ pl\k^U2-3jkU^-3U,\ 

+p 

+ pf[^^W,-3JkW,-3W,\} 


For the AH^ fields the following can be stated. 


AH^ = 




- V X AE^ = - 




-z,lV,+JkV,)] 


-f 


^\pAy,->-m)-p,(u,+ 

Spi 



(67) 

and simplifying yields the 

( 68 ) 

(69) 

(70) 

(71) 
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Substituting the partial derivatives of Un and Vn into AH^ and simplifying yields the 
following. [Ref 6 ] 


(AH^) +3F=) 

^ JCOHATt 

+ Aa(3^5 + 3 jkW^ - + 3 C /5 + 3 jku^ - k^u) } (72) 


Two additional methods (Method B-MaxwelTs curl equations applied twice to 
original vector potentials and Method C-solving AE^ and AH directly from vector 
potentials by an alternative method) for determining the cross-field terms are derived in 
Appendix-A. 
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m. CALCULATING FORWABD PROPAGATION FIELDS AND TESTING 

THEIR CONVERGENCE 


In this chapter, the AE and AH segments are used to generate a transfer matrix. 
The fields on a cylindrical surface at a distance of p 2 from the source can be determined 
by applying this transfer matrix to the equivalent currents generated by fields at a 
distance of pi from the source. The convergence of the three methods for determining 
AE and AH is examined. Choosing Method A, a parametric analysis of the convergence 
of fields at P 2 is conducted. Trends of the convergence of the calculated fields to the 
exact fields at P 2 are examined based upon differing input parameters in Method A. 


A. CALCULATION OF TRANSFER MATRIX FOR CALCULATING THE p2 
FIELDS 

Each EJ^, Ep , and/fj^ segment depends on several parameters, the most 
important of which is the distance R between the observation point on a cylindrical 
surface at a distance P 2 from the source and points on a cylindrical surface at pi. In order 
to estimate the fields at the i* point on p 2 , the contributions from all points on pi must be 
sununed. Figure 10 shows the physical relationship between the points on p 2 and the 
points on pi. 
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fe - ■£«]„ ; TO 

These vector multiplication equations can be condensed into a vector-matrix 
multiplication as shown below. 

^11 ^1/ •" ^li 

[^1 ••• ^n\pi • • * 

J’m • • • '^Ni • • • ^At 

TO 


= [£. ... (77) 

P\’-*P2 


If the point distributions on pi and p 2 are uniform, then all and all 

. If the point distributions are uniform and equivalent for both pi and P 2 , 
then Az„ = Az„ = Az. If either of the distributions on the pi or P 2 surfaces is not 

uniformly spaced, then the resultant transfer matrix T will have no commonly 
identifiable structure. If the point distributions are both uniformly spaced, but with 
tsz„ tiz „, then T will have some structure. When Az„ = Az^ , then the transfer 

matrix will be banded and have a Toeplitz structure [Figure 1 l][Ref 7], with t„ ^ t _„. If 
Az„ = Az^ and z„ = z„, then the transfer matrix will have a Toeplitz structure with 

Pi Pi Pi Pi ^ r' 



Figure 11. General Structure of a Toeplitz Marix 
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The bands in the transfer matrix are caused by the - z^^„) ^ terms in R, (AE^ 

and AH, and - z^^„) * terms in (AE)p. When Az’s are uniform and equal, but the 

actual z distributions are not equivalent (i.e. offset from each other), then the transfer 
matrix elements are translation invariant and depend only on m-n. When Az’s are 
uniform and the z distributions are aligned, then only the (AE)^term is translationally 

invariant depending on m-n. All the other terms are based on (z^^„ - z^„) \ so they are 

translationally invariant depending on Im-nl. Transfer matrices with the m-n 
dependence are full, nonsymmetric Toeplitz matrices, in Figure 11. Transfer 

matrices with the |m-ni dependence are full, symmetric, non-Hermetian Toeplitz 
matrices, t„ =t_„ in Figure 11. An interesting side note is that transfer matrices for 

Azi=CAz 2 or AziD=Az 2 , with C and D integer constants greater than one, are made by 
selecting specific rows or columns of a Toeplitz matrix based on the Az with the smallest 
size. 


The specific fields required for this thesis’ study of back propagation are the Ez 
and Hj, fields on the p 2 cylindrical surface. These total fields are found using a large 
(2Mx2N) transfer matrix consisting of four blocks. Each block is the (MxN) transfer 
matrix for a specific contribution to the total field from one of the two equivalent 
currents on the pi cylindrical surface. 


'■2Mx2N 


L^HE 




‘HHJ 


(79) 


( 80 ) 


In general the complex elements of The very large in magnitude compared to the rest 
of the matrix elements. Conversely, the elements of Teh are very small compared to the 
other elements. The relative magnitude of the Tee and Thh elements are roughly the 
same. This difference in relative magnitudes is due the small H-field at pi contributing to 
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the large E-field at pa for The, and, conversely, the large E-field at pj contributing to the 
small H-field at pa for Teh. 

Additionally, T will generally be rectangular since more points will be 
determined at pa than are measmed at pi. This is in preparation for over determined 
inverse systems. 


B. COMPARISON OF THREE METHODS FOR DETERMINING THE 
CALCULATED FIELDS AT pa 


The three different methods for determining the calculated fields at Pa are 
compared using the Least Squares Error (LSE) of each method’s field compared to the 
exact field. 

LS£, = j2^S^l%:iI;rt00 (81) 

V aSFexact) 

The LSE is a measure of a particular method’s total energy contained in the error with 
respect to the total energy contained in the exact field, expressed as a percentage. A 
comparison of the three methods LSE shown in Figure 12 reveals that their differences 
were indistinguishable. This indicates that no particular method is better than the others 
in the near field. Method A was chosen for all further investigations. For the fields 
calculated in Figure 12 and in all subsequent figures, the frequency of the current 
distribution on the dipole is /= 300 MHz. This frequency corresponds to a wavelength 
X = 1 meter. Thus, values of p and z also represent values relative to X. 
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Figure 12. LSEs for Three Methods on the Same Plot 











C. CUMULATIVE TOTAL ENERGY DENSITY OF FIELDS 


In order to get an idea about the spread of energy contained in the fields of a 
single dipole, the c um ulative total energy distribution for both Ez and H4, is determined in 
Figure 13 and Figure 14. For a dipole centered at the origin, the plots illustrate the 
percent of total field energy in the field for | z | >Zc, the z truncation. These plots provide 
a guideline for choosing a maximum z in order to keep a standard amount of field energy 
when doing calculations. Additionally, it is also a measure of the rapidness of change in 
percent of energy as Zc approaches zero. These two considerations provide a starting 
point for understanding the effects of changing parameters to the amount of error to be 
expected in the calculated fields at p 2 . 

D. PARAMETRIC ANALYSIS OF THE CALCULATED FIELDS USING 
METHOD A 

The following subsections outline trends in the accuracy of the calculated fields 
with respect to the exact fields on the cylindrical surface at P 2 . All calculations use 
uniformly distributed, equivalent point distributions for z^and z^ 2 > with the z^^ 

maximum/minimum equal to the z^^ maximum/minimum. The relative merit of various 

parameter changes can be reviewed using three techmques. First, the Real and Imaginary 
components of the exact and calculated fields are plotted together. Second, the Relative 
Error of the magnitude of the calculated field with respect to the magmtude of the exact 
field on a point to point basis is plotted. Last, the Least Squares Error expressed as a 
percentage is plotted. 
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Contour Plots of Cumlative Energy Distribution for Mag(E^) from Dipole at the Origin 



Figure 13. E Plot 
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p (m) 



Cutoff- z (m); zMax = 20 (m); Dipole Length = 0.1 (m); Nz = 20000 


Figure 14. H Plot 


31 




1. Varying the Truncation in Z 


To get an appreciation of the effect of varying the truncation of z at a fixed p 
value, Figure 15 shows a plot of the LSE in the calculated pz fields for a range of 
maximum z. From Figures 13 and 14 it is clear that as the cylindrical surface pi moves 
further from the source the value of the z cutoflF point must be increased to maintain a 
constant level of the field’s energy for pz calculations. Figure 13 shows that changes in 
the z truncation point have a significant impact on the amount of error in the pz fields 
when the magnitude of the maximum z is small. Then, beyond a certain point, increasing 
the truncation point in z does little to decrease the error in magnitude z pz fields. Figure 
15 confirms these two observations. For z cutoff less then two meters, the LSE is greater 
than one percent. 

2. Varying the Length of Az 


In building the transfer matrix and defining the fields on the two cylindrical 
surfaces, an estimation of the actual smooth continuous fields was made. The fields were 
estimated to be the same over a length of z, called Az. As Az decreases, higher spatial 
fi-equencies for fields can be resolved uniquely. In other words, the estimate of the 
continuous field made by discrete steps Az apart becomes more accurate as Az decreases. 
If there are no spatial frequencies higher than some kz, then decreasing Az when 

'2,7c 

Az « — will provide little improvement to the calculated field accuracy. Conversely, 
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decreasing Az when Az » — will dramatically improve the field accuracy. 

k.c 


These trends are easily seen for cases when the z truncation point is large enough 
to include most of the energy in the field. As seen in Figure 16, when coarse estimates of 
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the field are refined, large Az is decreased and the error in the calculated fields also 
decreases. 
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Least Squares Error (%) for Mag(E^) for Varying Max 2 and fj2; pi = 0.5 m; N^= 100 



Max z (m): Nz/meter = 32 


Least Squares Error (%) for Mag(H^ for Varying Max 2 and p2; pi = 0.5 m; 100 



Max 2 (m); Nz/meter=32 


Figure 15. Plots of Varying Maximum z 
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fi2(m) 


#2 points/meter: zMax= 1 


Figure 16. Plots of and for Vaiying Az 
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Two special cases exist, both of which are caused by rapid fluctuations and large 


estimates obtained when R becomes small, and — becomes large. The first case occurs 

R" 


throughout the region of interest when Ap=p 2 -pi becomes small. This causes the 
calculated fields to be very large and locally driven, even when these calculated fields are 
well behaved for a p 2 only a few more Ap lengths away. Decreasing the size of Az is 
required to compensate for this effect as seen in Figure 17. The second case occurs when 
both Pi and P 2 are small in general. This means that fields close to the source are being 
estimated. These fields have rapid variation in the z direction and have large magnitudes 
as well. Additionally, since pi is small and pz is small, pi-p 2 will also be small, thus 
exacerbating the error by including large calculated fields. Figure 18 shows that 
decreasing the size of Az again helps to counteract this effect. 


3. Varying the Number of Thin Ring Segments 


Varying the number of segments (N^) on the ring leads to changes in the segment 
size. The segment size is equal to PiA(J)i. 

PM = ^ (82) 

As increases, the length of the segment size decreases and a more accurate calculation 
for Un, Vn, Wn, and the fields is produced. As an initial estimate, the segment size should 
be small compared to of the testing frequency. A more conservative limit on segment 
size would be to make it less than or equal to the size of Az. Alternatively, the segment 
size can be limited to between a fifth and a tenth of the size of Ap=p 2 -pi. This limit 
provides insight into why segments should be at least this small if not smaller. As p 2 gets 
closer to pi, the points most influenced by the thin ring and by the ring segments become 
close to the center of the actual segment as shown in Figure 19. 
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Figure 18. Case II: pi Close to Zero, with Small Ap and Moderate Ap 
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Segment Area 


Large Ap 



Large Segment Area 



Large Segment Area 

Figure 19. Effect of Segment Size on Error 

If the overall size is not reduced, then that segment’s calculated value is not a 
good estimate for the field contributions it represents. Figure 20 illustrates the effects of 
small versus large N^, and large segment size versus small segment size, respectively. 
Clearly, for small the calculated fields are a poor estimate. Increasing N 4 , quickly 
improves estimates. For further calculations the segment length will be ~Az for small 
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Az. This gives N. = ?:EBl for small Az. For all calculations used in the plots, is set 
^ Az 

to or 100, whichever is the larger of the two. 

Az 

4. Overall Trends 


From a review of the previous figures, several other interesting features can be 
discussed. First, the calculated fields all have errors close to the z truncation points. 
This is clearly shown on the plots of Relative Error. This is caused by not including the 
influence on pi at points I z | > z cutoff By truncating fields greater than a certain z 
maximiim^ error is introduced. The error at the edges of the calculated field is due to the 
proximity of field points on pi that are not included when making an estimation of P 2 
fields. 

Another less obvious effect occurs when increasing p 2 from a fixed pi. After an 
initial dip in error, LSE plots show error slowly increasing again as P 2 increases. 
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(m): Nz/meter = 32 



Figure 20. Plots Showing Effect of Small Versus Large 
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IV. BACK PROPAGATION OF FIELDS 


In forward propagation of fields, the contributions of tWn rings of electric and 
magnetic equivalent currents determined from the fields at the pi cylindrical surface were 
summed to produce an estimate of the exact fields at the p 2 cylindrical surface. This 
process was described mathematically in the form of vector-matrix multiplication. 





jCALCULATED 


(83) 


Algebraically, to obtain the fields Ez and H 4 , at pi from knowledge of the calculated 
fields at P 2 , it appears that multiplying both sides of the equation by the inverse of the 
transfer matrix T will yield the correct answer. 

[£Ai,[7'][?’]-' = M.U7'r‘ (84) 




(85) 


Figure 21 shows the effect of T* on the calculated Ez and fields of p 2 . The reformed 
Ez and H* fields of pr tend to be very close to the exact ones. In fact, the difference is 
within machine accuracy in some areas. This is expected since the T matrix actually 
formed the calculated fields. 
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PoInMso Ret Error In Exact-Num Soln for Mag INV (p 1) Pointwieo Rel Error In Exact to Num Soln for Mag INV E^ (p 1) 


.5 Relative Error in Exact to Num Soln A for Mag INV E, fo 1) from p 2 Calc Total Reids and InvTT 

X10 ^ 



^ Relative Error in Exact to Num Soln A for Mag INV (p 1) from p 2 Calc Total Fields and InvTT 

10 r' 



Figure 21. Effect of T'^ on the Calculated Ez and Fields of P 2 . 
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A. EFFECT OF APPLYING INVERSE T MATRIX TO THE EXACT FIELDS 


ATp2 

When the calculated is replaced by the exact , a well behaved 

would yield newly calculated inverse fields close to the source fields at pi. If the 
difference between the calculated and exact fields at p2 is small, then appl5dng a well 
behaved T’ would yield a newly calculated field whose values would be close to the 
exact Pi field values. Mathematically, the following relationship holds true for a well 
behaved T* and 5 invt is small. 

+ ^LCAUZtoEXACT)^' Y 

Unfortunately, T'^ is not well behaved in this case. Inverse T is extremely sensitive to 
small changes in the fields at p2 used as its input. As can be seen from Figure 22 , even a 
small change of 5 a between the calculated and exact fields at p2 will cause large 
differences between the exact and inverted fields at pi. At first glance, it seems that the 
best way to get better results is to make the small change in 5 a even smaller. Figures 23 
and 24 show the pointwise relative error of the magnitude of the inverted fields compared 
to the exact fields at pi. Two cases are shown in each figure, one for a segment density 
of eight per meter and the other for a segment density of sixteen per meter. The LSE for 
the inverted pi fields and the calculated P2 fields is also shown to give a measure for the 
overall improvement increasing segment density provides. Also the maximum 5 a 
between calcualted and exact fields at p2 is provided. 

Part of the reason for the poor improvement in field esstimation lies with T and 
T'\ They are classically ill-conditioned matrices. Small changes at the input of an ill- 
conditioned matrix cause large changes at the output. One measure of this effect for a 
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matrix is called die condition number. When condition numbers are very large, matrices 
are considered ill-conditioned. Additionally, the transfer matrix becomes physically 
larger when either more terms are used to determine P 2 fields or when more p 2 field 
points are desired. Therefore, as the accuracy of the calculated fields improves, more 
points have been used and T becomes larger and larger. If T is an ill-conditioned matrix. 
Figure 25 shows that larger versions of T are even more ill-conditioned. The growing 
condition number is one reason that decreasing 5^ has only a small improvement in the 
accuracy of the calculated inverse fields. Therefore, in order to get a better, smaller 5 a , 
the size of T must be increased, but in increasing T, T*^ is increased and becomes even 
more ill-conditioned, countering the effects of a smaller 5 a. 
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(p 1) and INV (p 1) ^ ^ H^(p 1) and INV H^(p 1) 



z(p1) (m): Density N^, = 16; p 1 = 0.2 (m); p 2 = 0.3 (m); = 100 



z(p1) (m); Density = 16; p 1 = 0.2 (m); p 2 = 0.3 (m); = 100 


Figure 22. T'* Sensitivity. 
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Relative Error in Exact to Num Soln A for Mag INV E^ (p 1) from p 2 Exact Fields and InvTT 
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Figure 23. Effect of Decreasing 5 a on Error in Electric Fields 
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-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 


z(p1) (m): Density = 8; p 1 = 0.2 (m); p 2 = 0.3 (m); = 100 



z(p1) (m); Density = 16; p 1 = 0.2 (m); p 2 = 0.3 (m); = 100 


Figure 24. Effect of Decreasing 6^ on Error in Magnetic Fields 
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10 20 30 40 50 60 


#z points/meter; zMax = 2 tn; p2 = 0.3 iti; pi = 0.2 m; N^= 100 

Figure 25. Condition Numbers of the Transfer Matrix 
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B. USING AN OVERDETERMINED INVERSE T MATRIX 


In an attempt to improve the accuracy of the inverted fields, the inverse T matrix 
is constructed to be over determined. The fields at P 2 were determined at more points 
than the fields at pi. This leads to a “long and skinny” T matrix as indicated by the 
following equation. 

N>M 

Since T is no longer a square matrix, the inverse of T is not quite as straight forward. 
Instead of the standard A'^ of a nonsingular square matrix, the pseudoinverse, A"^, is 
calculated using the SVD technique shown in equation (88) where U and V are unitary 
matrices and S is a diagonal matrix. The pseudoinverse is calculated in equation 89 from 
the inverse ofS, S'*’. 

'^1 

A=USV' S= s 
_0 

Jl_ 

^1 

$2 

0 

Note, when A is 2Mx2N, with N > M, then is 2Nx2M and AA'*’=I( 2 Mx 2 M)- By using 
the SVD to find the pseudoinverse, the matrix T"^ will find the least squares estimate of 
the Pi Inverse Fields from the over determined inputs [Ref 9]. Figure 26 shows the LSE 
inverted fields at pi when increasing the number of P 2 segments by multiples of the 


A^=VS'tl' S* = 
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Least Squares Em>r (%); M|ig(INV H '-oast Squares Error (%); l^g(INV E 


number of pi segments. After the ratio of p 2 to pi segment densities increases beyond 
three, very little improvement in the inverted fields’ LSE is observed. In fact, the LSE 
degrades as the ratio increases from one to three. 






C. IMPROVING INTEGRATION ESTIMATE ON pi WHILE MAINTAINING 

A SMALLER T MATRIX 


The next attempt at improving the calculated inverse field is to maintain a smaller 
T matrix size while attempting to get some benefit from more integration accuracy. In 
order to do this, a new matrix D is defined such that it is a 2CN x 2M matrix created 
from a z distribution as shown in Figure 27. Then bands containing an odd number, C, of 
the matrix’s rows are added together and divided by the number of rows in the band, C. 
This new averaged integration row is then used as a row in T, a 2N x 2M matrix, where it 
is multiplied by the field point corresponding to the original D matrix band’s middle row. 


r 


C points of 
the D transfer 
matrix 



point used for Field Evaluation 
^^"^ntegration Contributions Averaged Together 


( 3 ) 




z point used for Field Evaluation 
^^'^ntegration Contributions Averted Together 


Figure 27. New Approximate Integration Row Used in T. 
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Thus, becomes ^\}rx 2 M by averaging bands of C rows, where C is and odd 

integer. Figure 28 shows the LSE of inverted fields at pi for an increasing ratio of pi 
segment density used to improve the integration accuracy to the P 2 segment density. As 
the ratio increases beyond three, the inverted field LSE changes only slightly. 


D. A CLOSER LOOK AT REAL AND IMAGINARY FIELD COMPONENTS 


Figures 29 and 30 show the Real and Imaginary E-field components for the exact 
Pi E-fields and the pi E-fields calculated using the matrix T inversion process for several 
Pi, P2 pairs. Figures 31 and 32 show the Real and Imaginary H-field components for pi 
H-fields. It is striking that in all cases, the Real Inverted E-fields and the Imaginary 
Inverted H-fields are essentially zero. 
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east Squ 


Least Squares Error (%) for Mag{lNV E^) for Varying Integration Nz Density 



Integration Density CONST; zWlax = 2 m; p1 = 0.2 m; p2 = 0.3 m; N^= 600 


Least Squares Error (%) for Mag(INV H^) for Varying Integration Nz Density 



Integration Density CONST; zMax = 2 m; p1 = 0.2 m; p2 = 0.3 m; N^= 600 


Figxire 28. Flat Weighted Integration—^Inverse Fields 

















2{p1) (m): Density = 32; p 1 = 0.2 (m); p 2 = 0.3 (m); = 100 


Figure 29. Back-Propagated and Exact Electric Field Components, Case 1 



2(p1) (m): Density = 32; p 1 = 0.5 (m); p2 = 0.6 (m); = 110 


Figure 30. Back-Propagated and Exact Electric Field Components, Case 2 


56 









Exact H^(p1) and Num Soln for INV H^(p1) 



z(p1){m): Density N^ = 32; p1 = 0.2(m); p2 = 0.3(m); N^= 100 


Figure 31. Back-Progagated and Exact Magnetic Field Components, Case 1 



z(p1) (m): Density = 32; p 1 = 0.5 (m); p 2 - 0.6 (m); - 110 

Figure 32. Back-Propagaged and Exact Magnetic Field Components, Case 2 
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The inverted fields seem to be composed of only the imaginary E-field and the 
Real H-field. To investigate the cause of this, the inverted fields at pi were propagated 
back to p 2 to validate the T matrix. As expected from a transform, inverse transform pair 
of operations, the inverted fields at pi recreated the calculated fields at P 2 . Thus, using 
the T matrix, there are at least two field distributions on pi that create the fields at pz. 
Upon further investigation, various combinations of the Real and Imaginary components 
of the exact fields on the pi cylindrical surface were propagated out to the pa cylindrical 
surface using the T matrix. Figures 33 through 36 show that the combination of the 
Imaginary E-field and Real H-field components of the pi field distribution recreated the 
desired fields at pa. The two field components that formed the calculated pa field 
distribution are also the only non-zero fields produced when T'* is applied to the exact pa 
fields. Figures 35 and 36 show that the combination of Real E-field and Imaginary H- 
field components of the pi field distributions seemed to produce zero fields at pa- This 
led to the discovery that the real pi E-field and the imaginary pi H-field are 
complementary solutions to Maxwell’s equations and are in the null space of the 
transform matrix T operator. [Ref 10] Conversely this implies that there is a particular 
solution to Maxwell’s equations that produces the fields at pa. [Ref. 10] 
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Figure 33. Exact and Particular Solution Component Calculated Magnetic Fields 


Exact and Particular Soln for E^(p2) 
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Figure 34. Exact and Particular Solution Component Calculated Electric Fields 
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Exact Ej,(p2) and Null Soln Ibr E^(p2) 



2(p2) (m); Density N^, = 32; p 1 = 0.2 (m); p 2 = 0.3 (m); = 100 


Figure 35. Exact and Null Solution Component Calculated Electric Fields 



2(p2) (m): Density = 32; p 1 = 0.2 (m); p 2 = 0.3 (m); = 100 


Figure 36. Exact and Null Solution Component Calculated Magnetic Fields 
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In order to better understand this. Maxwell’s equations are examined for a strictly 
real current and for E and H-fields broken into real and imaginary components. While 
J = + jj^ in general, the assumed dipole current distribution is strictly real, J = J^. 


Ea = Eaz + JEai 

(90) 

Ha- Ha\ + JHaz 

(91) 


The real and imaginary parts of Ea and Ha are indexed such that paired component 
solutions have identical indices. Substituting equations (90) and (91) into Maxwell’s 
equations illustrates the respective field components relations. 

V X ^^2 + JEai)= -jco^Ax + JHai) (92) 

V X {^Ai + jae^Ai + jEA))+JA (93) 

These two equations are then split into four by grouping real and imaginary terms. 


V X {jEa^= -ja)/n^A^ 

(94) 

Vx^^i)= jG)s{jEA\)+J A 

(95) 

V X (£^2)= -ja)fi{jHAz) 

(96) 

V X {jHaz) - j(OS^Az)+ 0 

(97) 


Thus £’.41 and Hai in equations (94) and (95) solve the particular solution for Maxwell’s 
equations with a real source current J a . Similarly, Eai and Haz in equations (96) and 
(97) solve the complimentary or source free solution for Maxwell’s equations. [Ref 10] 

The fields partitioned into particular and complimentary solutions of Maxwell’s 
equations can then be transformed into equivalent electric and magnetic current 
distributions by the surface equivalence theorem. The equivelent surface currents are 
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divided into those created by the particular fields and those created by the complimentary 
fields as shown in Figure 37. 



Figure 37. Equivalent Surface Currents for Partioned Field Solution 
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(98) 


JsA\ +JsA 2 — nX\H A 1 + jH A 2 

Msa\ + MsA2 = {^A2+jEA^h (99) 


From the Surface Equivalence Theorem, one can now find the influence of just the 
particular solution fields and propagated fields in Region 2 outside the equivalent 
surface. A similar analysis can be done for the complimentary solution fields. 


JsA and Msa generate 


JsAi = nx-HAi 
Msai = JEax X h 


generate 


E = 0 

H = 0 

Regionl 

E = Ea2 + JEai 

Region2 

H=HAl+jHA2 


E = Ea2 

H = jHA2 

Regionl 

E = Ea2 + JEax 

Regionl 

H=liAX+jHA2 



J SA2 = Ha2 

Msa 2— jEA2'xn 


generate 


E = —Ea 2 
H=-jHA2 
E = 0 
^ = 0 


Regionl 

Region2 


[Ref 10] 


These results match the results fi-om using the outward propagation transform 
matrix T and the particular and complimentary fields. The above equations suggest that 
using the transform matrix T on particular solution fields will give the complementary 
fields inside the surface. Put another way, use the T matrix on fields generated at P 2 to 
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get some informaiton about fields at pi. Figures 38 through 41 confirm this result. 
Figures 38 and 39 show the contributions to the fields inside pi from integrating the 
equivalent currents generated at pi by the particular solutions to Maxwell’s equations. 
Inside pi, the particular solution at pi (the Real H-field and the Imaginary E-field) create 
the null solutions inside pi (the Imaginary H-field and the Real E-field). Figures 40 and 
41 show the contributions to the fields inside pi fi'om integrating the equivalent currents 
generated at pj by the null solutions to Maxwell’s equations. Inside pi, the null solutions 
at Pi create the inverse of the null solutions inside pi. Thus, the total field from null and 
particular solutions will cancel inside pi. Unfortunately, this method of back propagating 
a portion of the field gives only complimentary field information. Therefore, the inverse 
T method is still needed for back propagation of the particular solution. 
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Exact E-(p2) and Particular Soln for E^(p2) 
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Figure 38. Exact and Integrated Particular Solution Electric Fields 


Exact H^(p2) and Particular Soln for H^(p2) 



Figure 39. Exact and Integrated Particular Solution Magnetic Fields 
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Exact Ej,(p2) and Null Soln for E^(p2) 



z(p2) (my, Density = 32; p 1 = 0.3 (m); p 2 = 0.2 (m); = 100 


Figure 40. Exact and Integrated Null Solution Electric Fields 



2((<2) (m): Density = 32; p 1 = 0.3 (m); p 2 = 0.2 (m); = 100 


Figure 41. Exact and Integrated Null Solution Magnetic Fields 
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V. CONCLUSION 


As the distance from a source increases, the ability to resolve rapid field 
variations close to the source degrades. Once in the far field, two sources can only be 
distinguished if they are greater than one half wavelength apart. If measurements are 
made on the surface of the object, the physical presence of the probe disrupts the fields it 
is attempting to sample. Moving a short distance away from the source into the near field 
provides a compromise between probe interference and resolution capability. The near 
field probe will no longer heavily disturb the fields and will be able to measure fields 
capable of providing better than one half wavelength localization of sources. 

The techniques used to analyze back propagation do not give a perfect 
reconstruction of fields closer to the source. In fact, some components of the calculated 
inverse fields are essentially zero and provide no useful information about the location of 
the source whatsoever. 

The reason that the inversion technique as used in this thesis does not give better 
answers is threefold. At the level of analyzing the transfer matrix structure, the matrices 
used are very ill-conditioned. In fact, as the matrices become more accurate—^more z 
range with finer integrations, they should become “perfectly” ill-conditioned. They 
transform non-zero vectors to zero exactly. The physical nature of the fields irmately 
creates ill-conditioned matrices. As more points and finer resolution are used to decrease 
error in the calculated fields at pz, the transfer matrix thus produced becomes more and 
more ill-conditioned. Inverting the ill-conditioned matrix creates another ill-conditioned 
transfer matrix for determining the source fields at pi. The matrix is so ill-conditioned 
that even a very small difference between the calculated and exact field at pz causes large 
error in the calculated inverse field. 

At an intermediate level, an actual radiator with a cylindrical surface at pi may 
produce evanescent fields as well as the fields identical to those from a real somce 
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current located inside pi. These evanescent fields would have a small but non-zero 
contribution to the fields measured at P 2 . For equivalent currents at pi, no actual 
evanescent fields exist. The equivalent currents are used to make numerical estimates of 
the fields at p 2 —the calculated P 2 fields. When inverting the exact fields at P 2 to estimate 
the fields at pi, the small differences between the calculated p 2 fields and the exact pz 
fields could be created by an evanescent field at pi. Thus, this technique has the 
potential to produce non-physical evanescent fields as well as the actual fields on pi. 

At the highest level, there are the cases of non-radiating fields [Ref 8 and 10]. 
The fields have a mathematical structure such that they can be large at pi and exactly 
zero at p 2 and at all points greater than pz. These fields are the homogeneous solution to 
the forward transfer matrix and can be considered to be in the null space of the forward 
transfer matrix. These fields can exist mathematically, but discretizing the problem may 
cause the inverse transfer matrix to reconstruct the actual along with some part of a non¬ 
radiating or evanescent field. 

Therefore, by its very nature the inverse source problem as handled in this 
technique is not only ill-conditioned, but is perfectly ill-conditioned, since null or very 
small fields at pz can exist due to large non-physical fields at pu. There will then be the 
actual inverse field plus any number of non-radiating and evanescent fields that can exist 
at Pi and still give the measured exact fields at pz. 

The technique was chosen for its ability to be expanded to any axisymmetric 
problem, not just cylindrical symmetries, in a manner that did not exclusively depend on 
a separable geometry and specialized Hankel functions. 

There are several areas for further study using this technique. The first is a study 
of modifications to the transfer function in order to eliminate its null space in the 
inversion process. The second is examining methods for decorrelating the terms in the T 
matrix to simplify inversioiL The third is wavelet analysis of improved transfer matrices, 
specifically focusing on a wavelet’s joint space/spatial scale localization abilities. The 
fourth is applying this technique to other axisymmetric fields and sources (such as a 
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spherical radiator, conical radiator, and other sources found on a surface of revolution). 
Finally, this technique can be used to find the actual surface currents on metal radiators 
and scatterers. 
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APPENDIX A. METHOD B-MAXWELL’S CURL EQUATIONS APPLIED 
TWICE TO ORIGINAL VECTOR POTENTIALS AND METHOD C-SOLVING 

ae7 and ah directly from vector potentials by an 

ALTERNATIVE METHOD 


A. Method B: Maxwell’s Curl Equations Applied Twice to Original Vector 
Potentials 

In the previous section, AH^ and Af"^ were determined from AE^ and AH^ 
respectively. In this method, AHj? and AE^ are found directly without the intermediate 
step of completely determining and AH^. Since 


=-JLvxiJ = — 


J(oe, 


then = ■ 


1 


jo>s^ 


Vx 


JG>S, 

1 


Vx 

' 1^ 
X 

> 

1 

1_ 


JJ 


Vx 


il 


IR 




R 


-ds' 


( 100 ) 


( 101 ) 


The curl operations can both be brought inside the integral to give the following. 

E^=—i—f Vx Vx J— (102) 

JQ)SMs I L V ^ /J. 

Breaking the surface of the integration into many thin rings gives the following. 


AE, 


^ A^i f 

jasMh 


TT 





f 

Vx 

Vx 

jf 


J 

-7t 



[ R ). 





(103) 


In section n.B.2 ‘The Equivalent Electric Current Case and the Vector A Potential,” it 
was found the following was found. 
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(104) 


VxJJ - 

^ R 


(-1 -jk 
R^ 



Substituting this term yields the following. 


Vx 







X 


rj_ 

W 




This can be further expanded by using the 



functions as follows. 


Vx 


Vx7,z 



Ll + M 

^R^ R^ 



Solving for x [r x z| and [v x [r x z| and substituting Un, V^, and Wn 
gives the final result. 

(AEj = ^ fe ^ 

^ ja €, 4^ 

+p,^%-3JkV,-3V,^ 

{&E,l = -3JkU,-3U,) 

jco e, 4;r 

+ p,P2(6V,+6jW,-2k%)\ 


(105) 

(106) 

functions then 

(107) 

(108) 
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Note that while (AE^) is the same for Method A and Method B, (AE^)^ is somewhat 


different Similarly, the following is true. 


= Vx VxF 


( 109 ) 


Hp=- 


1 „ 1 ^ 


Vx Vx —\\M,- - ds’ 


( 110 ) 


Hp “ 


jCOfl.AK 


JJ Vx VxA< 


( 111 ) 


Once again, breaking the cylindrical surface at pi into thin rings will give AH. 


e--'“'(vx(Fx|^cos(iJi -^sin^Ji)))) 


( 112 ) 


(AH^), = k + >^2 + fe -- VkV, -3V,) 

+ + 3 jkU^ - + 3 ^ 5+3 jkW^ - k‘^W^\ ( 113 ) 


This is the same (AH^V as determined in Method A. 


75 



B. Method C: Solving AE^ and AH Directly From Vector Potentials by an 
Alternative Method. 

In the previous section AE^ and AH^? were solved directly jfrom Maxwell’s curl 

equations E. = ——V x and Hp = ——V y-Ep . An alternative method is to 
jG>s^ jeoM. 

substitute = —V x A and Ep = —V x F and use Lorentz condition. 


E.=-ja}A —^ v ( v * a ) 


Hp=-ja)F —^v(v*f) 

COJUS 




4;r " ‘ R 


- ds'- ^ 


G>sAn 


JJ 


-JkR\ 


V*J^Z- 


R 




W 




-JkR 


J 


. -ds’- 

An R (ojuAn 


III 


( — 
V V*M- 


R 


JA 


W 


Expand the second terms of E^ and Hp 


f 

V*JJ- 

V ^ J 


= V 


/ 


-jkR\ 


J,z*V— 

V ^ J 


-jkR 


R 


-(V./,z) 


v*yj 

= v 


e-^'^(jJ»R) 

^ J 


\R^ R^ J 

\ * / 


(114) 

(115) 

(116) 

(117) 

(118) 

(119) 
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+ + (123) 

Break the cylindrical surface into thin rings to calculate Hp and and substitute Un, 
Vn, and Wn functions. 

(AE^)^ is the same as in the previous two methods. 

+—lz^-z,i^,*VkU,-eu^-U,-JkU ^]} (124) 
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jAn 




+—Lp 2 ( 3 ir,+ 3 jm^ -Ft/3-3^5 - 3 jm^+Jt%)-V3 - jkv ^]} (125) 

cos. 


Note that while (AE^)^ is the same, (AE^)^ and (A//^)^ are somewhat different in form 
when compared to the previous two methods. 
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APPENDIX B. COMPUTER SOURCE CODES DEVELOPED 


% TestSa.m 

% By M.A, Morgan 13 June 97 

% Mod“8a on 30 May 98 substitutes ifft for fft and fft for ifft 
% to correct kz polarity. 

% Propagates exact H_phi, E__z and E_rho fields from a z-axis 
% array of specified dipoles between rhol and rho2 using 
% cylindrical wave function transfer functions of complex FFT’s. 

% Fields are plotted on a user specified 2-D grid in rho and z. 

% An N-element uniform array of identical dipole elements is the 
% source. A sequential phase shift on element input currents 
% provides beam steering of the specified array. 

% Dipole element currents have sinusoidal form of I(z)=Im sin(h"lz|) 
% which has exact near-field integration results found in 
% Jordan & Balmain pp 333-337 or Elliott pp 329-332. 

clear all; eta=120*pi; 

f=300;% f=input(’Enter frequency (MHz): '); 

L=input (* Enter dipole element length L (meters): M; 
h=L/2; 

wl=300/f; k=2*pi/wl; k2=k*k; kh=k*h; Im=l; ckh=cos(kh); 

N=input(’Enter number of dipoles in array: ’); 
if N > 1, 

d=input{*Enter array spacing (meters): ’); 
phi=input('Enter far-field pointing angle (deg): '); 
zd=-(N-1)*d/2:d:(N-1)*d/2; % Element positions 

alpha=-k*zd*cos(pi*phi/180); % Element phasing 
I=Im*exp(j*alpha); % Input currents 

else, 

zd=0; I=Im; 

end 

rhol=input{'Enter rhol (m) for field calculation grid: '); 
rho2=input('Enter rho2 (m) for field calculation grid: ')/ 
nr=input('Enter Number of rho-Points: '); 
rho=linspace(rhol,rho2,nr); % Row array 

zl=input('Enter zl (m) for field calculation grid: '); 
z2=input{'Enter z2 (m) for field calculation grid: '); 
nz=input{'Enter Number of z-Points (64, 128, 256, ... etc): ')/ 
dz=(z2-zl)/(nz-1); zee={zl:dz:z2)'; % Column array 

% Wavenumber spectra 

dkz=2*pi/(nz*dz); kz2=nz*dkz/2; kzl=-kz2+dkz; 
kz=kzl:dkz:kz2; kzp=0:dkz:kz2; nkzp=length(kzp); 
nsup=find(kzp <= k); nsub=find(kzp > k); 
krsup=sqrt(k2-kzp(nsup).*kzp(nsup)); 
if '^isempty (nsub), 

krsub=sqrt(kzp(nsub).*kzp(nsub)-k2); end 
% Initializing Arrays 
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FEl=zeros(nkzp, 1); FH1=FE1; Hp=zeros(nz,nr); Ez=Hp; Er=Hp; 
[Rho, Zee] =itieshgrid(rho, zee); Rho2=Rho.*Rho; 


for n=l:N % Array element index 

ZnO=Zee-zd(n); Znl=ZnO-h; Zn2=ZnO+h; 

R0=sqrt(Rho2 + Zn0.''2); Rl=sqrt (Rho2 + Znl.''2); R2=sqrt (Rho2 + 

Zn2.^2); 

E0=2*exp(-j*k*R0); El=exp(-j*k*Rl); E2=exp(-j*k*R2); 

Hp(:,:)=Hp{:,:)+j*I(n)*(El + E2 - ckh*E0)./(4*pi*Rho); 

Ez(:,:)=Ez{:, :)-j*30*I(n)*{El./Rl + E2./R2 - ckh*E0./R0); 

Er(:,:)=Er(:,:)+j*30*I(n)*... 

(Znl.*El./Rl + Zn2.*E2./R2 - ckh*ZnO.*EO./RO)./Rho; 

end 

svfl=input(’Enter 1 to Save E and H Field Files: ’); 
if '^'isempty(svfl) ^ if svfl == 1, 

hdr=[f rhol rho2 nr zl z2 nz]; 
save hdr.dat hdr /ascii 
dum=real(Hp); save hpr.dat dim -ascii 
dum=imag(Hp); save hpi.dat dum -ascii 
duiri=real (Ez); save ezr.dat dum -ascii 
dum=imag(Ez); save ezi.dat dum -ascii 
dum=real(Er); save err.dat dum -ascii 
dxm=imag(Er) ; save eri.dat dum -ascii 
end; end 

F^l= 2 eros(nz,nr); Ppl=zeros(nz,1); Pp2=Ppl; Pp3=Ppl; 

% FFT in z of each fixed rho column 
% using ifft in z to correct polarity of k_z in FFT 
Hpf=nz*ifft(Hp) ; E 2 f=nz*ifft(Ez); Erf=nz*ifft(Er); 

Pr=-pi*Rho.*real(Ez.*conj(Hp)); % Time-average rho-directed power/dz 

Pz=+pi*Rho.*real{Er.*conj(Hp)); % Time-average z-directed power/dz 

Prf=-Rho.*real(Ezf.*conj(Hpf))/2; % Power Spectral Densities (PSD) 
Pzf=+Rho.*real(Erf.*conj(Hpf))/2; 


while If 

disp(’Program 


disp (’ 
disp (’ 
disp(’ 
disp (’ 
disp (’ 
disp (’ 
disp(’ 
disp(’ 
disp (’ 
disp (’ 
disp (’ 


“> 

==> 

==> 


Options: ’) 

Stop Program’) 

H_phi 3-D Plot’) 

E_z 3-D Plot’) 

E_rho 3-D Plot’) 

P_rho 3-D Plot’) 

P_z 3-D Plot’) 

P-vector Plot’) 

Select Rhol (Before Option 8)’) 

Fixed Propagation from Rhol to Rho2’) 
Propagation Errors from Rhol through Rho2’) 


nplot=input(’Select Numbered Program Option: ’); 
if isempty(nplot)I nplot > 9 | nplot <=0, break; 


end 


if nplot==l, 

clf reset; surfl(Zee,Rho,abs(Hp),[0 0]) 
if N > 1, 

title([’|H_{phi}I for N=’,int2str(N),’ Array; d=’,num2str(d),.•. 
’ m; L=’,num2str(L),’ m; phi=’,num2str(phi),’deg; f=’,... 
num2str(f),’ MHz’]) 
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€ 1S @ 

title(['lH {phi}l for L=' , niim2str (L), ’ m Dipole at f=’,num2str (f), 
’ MHz']) 

snd 

xlabel(*2 (m) M ; ylabel('rho (m) ’); zlabeK’|H_{phi} I (A/m) *) 

view{200/30); v=axis; axis([zl z2 rhol rho2 v(5) v(6)3); figure(1) 

hepy=input(*Enter 1 for Hard Copy: ’); 

if -isempty(hcpy), if hepy == 1, print; end; end 

shading interp; view(90,-90); colorbar; figure(1) 

hepy=input(’Enter 1 for Hard Copy: *); 

if is empty (hepy) / if hepy == 1/ print; end; end 

elf reset; surf(Zee,Rho,loglO(abs(Hp))) 

if N > 1, , . ^ ^ 

title([’|H {phi}l for N=’,int2str(N),’ Array; d=’,num2str(d), ... 

’ m; L=’,num2str(L),’ m; phi=’,num2str(phi),’deg; f=’,..- 
num2str(f),’ MHz’]) 

e i se 

title(['|H_{phi}I for L=' ,num2str(L),' m Dipole at f=',nuin2str(f) / 
' MHz']) 

end 

xlabel('z (m) '); ylabel('rho (m)') 
v=axis; view(90,-90); shading interp 

axisdzl z2 rhol rho2 v(5) v(6)]); colorbar; figured) 
hcpy=input('Enter 1 for Hard Copy: '); 
if -isempty(hepy), if hepy == 1/ print; end; end 


% Manual fftshift for 1-D fft array 
Fpl (l:nkzp-2, :) =abs (Hpf (nlczp+1 :nz, :)); 

Ep>l (nkzp-l:nz, :)=abs (Hpf (Irnkzp, :)) ; 
elf reset; surfl (kz, rho, El)l', [0 0]) 
if N > 1, 

title(['H_{phi} Spectrum for N=',int2str(N), Array; 
d=', niun2str (d), ... 

' m; L=',nxm2str(L),' m; phi=’,num2str(phi),'deg; 
nvm2str(f),' MHz']) 


JL , • • • 


e i se 

title([’H_{phi} Speetrum for 
num2str(f),’ MHz’]) 


L=’,num2str(L),’ m Dipole at f=’. 


xlabel (’ k z (rad/m) ’); ylabel (’ rho (m) ’); zlabel (’ 1 FFT [H__^{phi} ] I ) 

view(200,30); v=axis; axis([kzl k 22 rhol rho2 v(5) v(6)]); figured) 

hepy=input(’Enter 1 for Hard Copy: ’); 

if '^isempty (hepy), if hepy == 1, print; end; end 

shading interp; view (90,-90); eolorbar; figure(l) 

hepy= input (’Enter 1 for Hard Copy: M; 

if '-isempty (hepy), if hepy == 1/ print; end; end 


elf reset; surf(kz,rho,loglO(Fpl’+•01)) 
if N > 1, 

title([’H_{phi} Speetrum for N=’,int2str(N),' Array; 
d=’,num2str(d) , ... 

* m; L=’,num2str (L), ’ m; phi=’,num2str (phi),’deg; f= ,... 
num2str(f),’ MHz’]) 

else, 

title([’H {phi} Speetrum for L=’,num2str(L),’ m Dipole at f=’,-.. 
num2str{f),’ MHz’]) 
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end 

xlabel('k_z (rad/m)'); ylabel('rho (m)•) 
v=axis; view(90,-90); shading interp 

axis([kzl kz2 rhol rho2 v(5) v(6)]); colorbar; figure(1) 
hcpy=input('Enter 1 for Hard Copy: '); 
if ~isempty(hcpy), if hcpy == 1, print; end; end 
end 

if nplot==2, 

clf reset; surfl(Zee,Rho,abs(Ez),[0 0]) 
if N > 1, 

title ([' |E_z I for N=', int2str (N), ' Array; d=’,nuia2str (d), ' la; L=’,— 
nuirL2str (L),' m; phi=',num2str (phi),'deg; f=',nuin2str (f), ' MHz']) 

else, 

title(['|E_z I for L=',num2str(L),' m Dipole at f=',num2str(f),' 

MHz']) 
end 

xlabel ('z (iti) '); ylabel ('rho (m) '); zlabel (' |E_z | (V/m) ') 

view(200,30) ; v=axis; axis([zl z2 rhol rho2 v(5) v(6)]); figured) 

hcpy=input('Enter 1 for Hard Copy: '); 

if ~isempty(hcpy), if hcpy == 1, print; end; end 

shading interp; view(90,-90); colorbar; figure(l) 

hcpy=input{'Enter 1 for Hard Copy: '); 

if ~isempty(hcpy), if hcpy == 1, print; end; end 

clf reset; surf(Zee,Rho,loglO(abs(Ez))) 
if N > 1, 

title (['|E_z 1 for N=', int2str (N), ' Array; d=',nuiu2str (d) , . .. 

' la; L=',num2str (L), ' m; phi=', num2str (phi),'deg; f=',— 
nioinZstr (f), ' MHz']) 

else, 

title ([' |E_z I for L=',num2str (L), ' m Dipole at f=',nuin2str (f),- 

' MHz']) 

end 

xlabel('z (m)'); ylabel('rho (m)') 
v=axis; view(90,-90); shading interp 

axis([zl z2 rhol rho2 v(5) v(6)]); colorbar; figured) 
hcpy=input('Enter 1 for Hard Copy: '); 
if ~isempty(hcpy), if hcpy == 1, print; end; end 

% Manual fftshift for 1-D fft array 
Fpl(1:nkzp-2,:)=abs(Ezf(nkzp+1:nz, :)); 

Fpl(nkzp-1:nz,:)=abs(Ezf(1:nkzp,:)); 
clf reset; surfl(kz,rho,Fpl',[0 0]) 
if N > 1, 

title (['E_z Spectrim for N=', int2str (N), ' Array; d=',nuia2str (d),— 

' m; L=',num2str (L), ' m; phi=',num2str(phi),'deg; f=',... 
num2str(f),' MHz']) 

else, 

title(['E_z Spectrum for L=',num2str(L),' m Dipole at f=',... 
num2str(f),' MHz']) 

end 

xlabel('k_z (rad/m)'); ylabel('rho (m)'); zlabel('|FFT [E_z]1') 
view(200,30); v=axis; axis([kzl kz2 rhol rho2 v(5) v(6)]); figured) 
hcpy=input{'Enter 1 for Hard Copy: '); 
if ~iseiapty (hcpy), if hcpy == 1, print; end; end 
shading interp; view(90,-90); colorbar; figured) 
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hcpy=input('Enter 1 for Hard Copy: 

if ~isempty(hcpy)/ if hcpy == 1, print; end; end 

clf reset; surf(kz,rho,loglO(Fpl' + . 01)) 
if N > 1, 

title(t'E_z Spectrum for N=',int2str(N), ' Array; d=',num2str(d), . 
' m; L=',num2str(L),' m; phi=',num2str(phi),'deg; f=',... 
num2str(f),' MHz']) 

else, 

title(['E_z Spectrum for L=',num2str(L),' m Dipole at f=',... 
n\mi2str(f), ' MHz']) 

end 

xlabel('k_z (rad/m)'); ylabel(|rho (m) ') 
v=axis; view(90,-90); shading interp 

axis{[kzl kz2 rhol rho2 v(5) v(6)]); colorbar; figure(l) 
hcpy=input('Enter 1 for Hard Copy: '); 
if ~isempty(hcpy), if hcpy == 1, print; end; end 
end 

if nplot==3, 

clf reset; surf1(Zee,Rho,abs(Er), [0 0]) 
if N > 1, 

title(['|E_{rho}I for N=',int2str (N), ' Array; d=', num2str(d),... 

' m; L=',num2str(L),' m; phi=',num2str(phi),... 

'deg; f=',num2str(f),' MHz']); 

else, 

title(['|E_{rho}I for L=',num2str(L), ' m Dipole at f=',num2str(f) 
' MHz']) 

end 

xlabel Cz (m) '); ylabel ('rho (m) '); zlabelC |E_{rho} | (V/m) ') 

view(200,30) ; v=axis; axis([zl z2 rhol rho2 v(5) v(6)]); figured) 

hcpy=input ('Enter 1 for Hard Copy: M; 

if ~isempty (hcpy), if hcpy == 1, print; end; end 

shading interp; view(90,-90); colorbar; figure(1) 

hepy=input('Enter 1 for Hard Copy: '); 

if ~isempty(hcpy), if hcpy == 1, print; end; end 

clf reset; surf(Zee,Rho,loglO(abs(Er))) 
if N > 1, 

title(['|E_{rho}I for N=',int2str(N),' Array; d=',num2str (d), ... 

' m; L=',niua2str(L),' m; phi=',num2str(phi),'deg; f=',... 
num2str(f),' MHz']) 

else, 

title(['|E_{rho]1 for L=',num2str(L),' m Dipole at f=',num2str(f) 
' MHz' ]) 

end 

xlabel('z (m) '); ylabel('rho (m)') 
v=axis; view(90,-90); shading interp 

axis([zl z2 rhol rho2 v(5) v(6)]); colorbar; figure(l) 
hcpy=input('Enter 1 for Hard Copy: '); 
if ~isempty(hcpy), if hcpy == 1, print; end; end 

% Manual fftshift for 1-D fft array 
Fpl(l:nkzp-2,:)=abs(Erf(nkzp+l:nz,:)); 

Fpl{nkzp-l:nz,:)=abs(Erf(l:nkzp,:)); 
clf reset; surfl(kz,rho,Fpl',(0 0]) 
if N > 1, 
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title([*E_{rho} Spectrum for N=',int2str(N),’ Array; 
d=’,num2str(d)__ 

* m; L=’,num2str(L),’ m; phi=',num2str(phi),'deg; 
num2str(f),' MHz']) 

else, 

title ([’E__{rho} Spectrum for L=',num2str (L) , ' rci Dipole at f=',... 
nuin2str (f), ' MHz']) 

end 

xlabel('k_z (rad/m)'); ylabel('rho (m) ') ; zlabel('|FFT [E_{rho}]|’) 

view(200,30); v=axis; axis([k 2 l kz2 rhol rho2 v(5) v(6)]); figure (1) 

hcpy=input('Enter 1 for Hard Copy: '); 

if -isempty(hcpy), if hcpy — 1, print; end; end 

shading interp; view(90,-90) ; colorbar; figured) 

hcpy=input('Enter 1 for Hard Copy: '); 

if -isempty(hcpy), if hcpy == 1, print; end; end 

clf reset; surf(kz,rho,logl0(Fpl'+.Ol)) 
if N > 1, 

title(['E_{rho} Spectrum for N=',int2str(N),' Array; 
d= ', nuin2str (d),. •. 

’ m; L=’,nm2str (L), ’ m; phi=’,num2str (phi),'deg; f=',... 
num2str(f),' MHz’]) 

else, 

title (['E_{rho} Spectriim for L=',num2str (L), ' m Dipole at f 
nu!ii2str (f) , ’ MHz']) 

end 

xlabel (' k __2 (rad/m) ') ; ylabel (’ rho (m) ’) 
v=axis; view(90,-90); shading interp 

axis([kzl k22 rhol rho2 v(5) v(6)]); colorbar; figure(1) 
hepy=input('Enter 1 for Hard Copy: '); 
if --isempty (hcpy), if hcpy == 1, print; end; end 
end 

if nplot~4, 

clf reset; surf1(Zee,Rho,Pr,[0 0]) 
if N > 1, 

title(['P_{rho} for N=',int2str(N), ’ Array; d=’,num2str(d),... 

' m; L=’,num2str (L),' m; phi=',num2str(phi),... 

'deg; f=',num2str(f),’ MHz’]); 

else, 

title([’P_{rho} for L=’,num2str(L), ' m Dipole at f=',num2str(f),.. 
’ MHz’]) 

end 

xlabel('z (m) ’); ylabel('rho (m) '); zlabel('P_{rho} (W/m)') 

view(200,30); v=axis; axis([zl z2 rhol rho2 v(5) v(6)]); figure(l) 

hcpy=input{'Enter 1 for Hard Copy: '); 

if -isempty(hcpy), if hcpy == 1, print; end; end 

shading interp; view(90,-90); colorbar; figure(l) 

hcpy=input('Enter 1 for Hard Copy: '); 

if -'isempty (hcpy), if hcpy == 1, print; end; end 

clf reset; surf(Zee,Rho,loglO(abs(Pr))) 
if N > 1, 

title(['P_{rho} for N=’,int2str(N),' Array; d=',num2str(d) , .. • 

’ m; L=',num2str(L),' m; phi=’,num2str(phi),'deg; f=’,..* 
num2str(f),’ MHz']) 

else. 
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title(['P_{rho} for L=',num2str(L),' m Dipole at f=*,num2str(f),... 
’ MHz’]) 

end 

xlabel('z (m) ') / ylabel{'rho (m)') 
v=axis; view(90,-90); shading interp 

axis([zl z2 rhol rho2 v(5) v(6)]); colorbar; figured) 
hcpy=input('Enter 1 for Hard Copy: '); 
if ~isempty(hcpy), if hcpy == 1, print; end; end 

% Manual fftshift for 1-D fft array 

Fpl(l:nkzp-2,:)=Prf(nkzp+l:nz, :); Fpl(nkzp-l:nz,:)=Prf{1:nkzp,:); 
clf reset; surfl(kz,rho,Fpl',[0 0]) 
if N > 1, 

title(['P_r PSD for N=',int2str(N),' Array; d=',num2str(d),... 

' m; L=',num2str(L),' m; phi=’,num2str(phi),'deg; f=',... 
num2str(f),' MHz']) 

else, 

title(['P_r PSD for L=',num2str(L),' m Dipole at f=',... 
num2str(f),' MHz']) 

end 

xlabel('k_z (rad/m)'); ylabel{'rho (m) ’); zlabel('P_r PSD') 

view(200,30); v=axis; axis([kzl kz2 rhol rho2 v(5) v(6)]); figure(1) 

hcpy=input('Enter 1 for Hard Copy: '); 

if ~iseiapty{hcpy), if hcpy == 1, print; end; end 

shading interp; view(90,-90); colorbar; figured) 

hcpy=input('Enter 1 for Hard Copy: '); 

if -isempty(hcpy), if hcpy == 1, print; end; end 

clf reset; surf(kz,rho,loglO(abs(Fpl)'+.01)) 
if N > 1, 

title(['P_r PSD for N=',int2str(N),' Array; d=',nuia2str(d),... 

' m; L=',num2str(L),' m; phi=',num2str(phi),'deg; f=',... 
num2str(f),' MHz']) 

else, 

title(t'P_r PSD for L=',num2str(L),' m Dipole at f=',... 
nuiti2str (f), ' MHz']) 

end 

xlabel('k_z (rad/m)'); ylabel('rho (m) ') 
v=axis; view(90,-90); shading interp 

axis([kzl kz2 rhol rho2 v(5) v(6)]); colorbar; figured) 
hcpy=input('Enter 1 for Hard Copy: '); 
if -isempty(hcpy), if hcpy == 1, print; end; end 
end 

if nplot==5, 

clf reset; surfl(Zee,Rho,Pz,[0 0]) 
if N > 1, 

title(['P_z for N=’,int2str(N),' Array; d=',num2str(d),... 

' m; L=',num2str(L),' m; phi=',num2str(phi),.., 

'deg; f=',num2str(f),' MHz']); 

else, 

title(['P_z for L=',num2str (L), ' m Dipole at f=', n\3iti2str (f), ... 

’ MHz']) 

end 

xlabel('z (m)'); ylabel('rho (m) ') ; zlabel('P_z (W/m)') 
view(200,30); v=axis; axis([zl z2 rhol rho2 v(5) v{6)]); figured) 
hcpy=input('Enter 1 for Hard Copy: '); 
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if ~isempty(hcpy), if hcpy == 1, print; end; end 
shading interp; view(90,-90); colorbar; figure(l) 
hepy=input('Enter 1 for Hard Copy: '); 
if ~isempty(hcpy), if hcpy == 1, print; end; end 

clf reset; surf(Zee,Rho,loglO(abs(Pz))) 
if N > 1, 

title (['P_z for N=', int2str (N), ' Array; d=',nuia2str (d) , ... 

' m; L=’,nimi2str (L), ' m; phi=',num2str (phi),'deg; f=',— 
num2str(f),' MHz']) 

else, 

title (['P_z for L=',nuiti2str (L), ' m Dipole at f=',num2str (f), ... 

' MHz' ]) 

end 

xlabel('z (m) *); ylabel('rho (m) ') 
v=axis; view(90,-90); shading interp 

axis([zl z2 rhol rho2 v(5) v(6)]); colorbar; figure(1) 
hcpy=input('Enter 1 for Hard Copy: '); 
if ~isempty(hcpy), if hcpy == 1, print; end; end 

% Manual fftshift for 1-D fft array 

Fpl (l:nlczp-2, : )=Pzf (nlczp+l:nz, :); E^)l (n]czp-l:nz, : )=Pzf (l:nlczp, :); 
clf reset; surf 1 (Icz, rho, E^l', [0 0]) 
if N > 1, 

title(['P_z PSD for N=',int2str(N),' Array; d=',num2str(d),... 

' m; L=',n:xiTi2str (L), ' m; phi=',nuia2str (phi),'deg; f=',... 
n\am2str (f), ' MHz']) 

else, 

title(['P_z PSD for L=',num2str(L),' m Dipole at f=',... 
nuiti2str (f), ' MHz']) 

end 

xlabel ('k_z (rad/ia)'); ylabel ('rho (m) '); zlabel('P_z PSD') 

view(200,30); v=axis; axis([lczl kz2 rhol rho2 v(5) v(6)]); figure (1) 

hcpy=input('Enter 1 for Hard Copy: '); 

if -isempty(hcpy), if hcpy == 1, print; end; end 

shading interp; view(90,-90); colorbar; figure(l) 

hcpy=input('Enter 1 for Hard Copy: '); 

if -isempty(hcpy), if hcpy == 1, print; end; end 

clf reset; surf (lcz,rho,logl0 (abs (Fpl)' + .01)) 
if N > 1, 

title(['P_z PSD for N=',int2str(N),' Array; d=',nuin2str(d),... 

' m; L=',num2str(L),' m; phi=',num2str(phi),'deg; f=',— 
num2str(f),' MHz']) 

else, 

title(['P_z PSD for L=',num2str(L),' m Dipole at f 
nuia2str (f), ' MHz']) 

end 

xlabel (' lc_z (rad/m) '); ylabel (' rho (m) ') 
v=axis; view(90,-90); shading interp 

axis([kzl kz2 rhol rho2 v(5) v(6)]); colorbar; figured) 
hepy=input('Enter 1 for Hard Copy: '); 
if -isempty(hcpy), if hcpy == 1, print; end; end 
end 

if nplot==6, 

npwr=input('Enter inverse scaling factor for quiver plot: ') 
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Pnnax=max (max (abs (Pr))); P 2 max=max (max (abs (Pz))); 

Pmax=min (Prmax, Pzmax); 

PrO=sign(Pr).*(abs(Pr)/Pmax).^(1/npwr); 

PzO=sign(Pz) .* (abs (Pz)/Pmax) .''(1/npwr); 
clf reset; quiver(Rho,Zee,PrO,PzO) 

title(['P~vector for N=’,int2str(N),' Array; d='/numZstr(d),.•. 

' m; L=',num2str(L),' m; phi=',num2str(phi),'deg; 
num2str(f),' MHz']) 

01S 0 

titl 0 ([ *P“V 0 Ctor for L=’,nuiti2str (L) / ’ m Dipolo at f=’,niini2str (f), ... 
» MHz»]) 

©nd 

xlaboK’rho (m) ’) ; ylabol (' z (m) ’) 

axis([rhol rho2, zl z2]); axis squar©; figur 0 (l) 

hcpy=input{'Entor 1 for Hard Copy: ]); 

if '^isoiapty (hcpy) / if hcpy ~ 1/ print; ©nd; end 

©nd 

if nplot“7, 

rl==input (’Enter Rhol radius to propagate from: ’); 

% Selecting closest grid radius 

[rm, ml]=min(abs(rho-rl)); rpl=rho(ml); 

% Forming denominator of radial transfer functions 
FEl=zeros(nkzp,1); FH1=FE1; 

[H, DH]=dhan)cel(0,krsup*rpl) ; FEl(nsup)=H; FHl (nsup)=DH; 
if ~isempty(nsub), , , , 

[K, DK]=dmodbes(0,krsub*rpl); FEl(nsub)=K; FHl(nsub)=DK; end 

% Energy Normalization for Optimal Estimators 
EEl=FEl'*FEl/nkzp; EH1=FH1'*FHl/nkzp; 

% Plotting Spatial and Spectral Fields at rho=rpl 
clf reset 

subplot(2,1,1), , , 

plot (zee, cibs (Ez (: ,ml)),'g', zee, abs (eta Hp(t,ml)), ;r ,... 

zee,abs(Er(:, ml)),'—c') 
if N > 1, 

title(['Fields and Spectra; Array N—',int2str(N), ; 
d=',num2str (d),— , 

'm; L=',num2str(L),'m; phi=',num2str(phi),'deg; f=',... 
num2str(f),'MHz; rho=',num2str(rpl),'m']) 

©Is© 

title(['Fields and Spectra for L=',num2str(L),'m Dipole at f=',... 
num2str(f),'MHz; rho=',num2str(rpl),'m']) 

end 

xlabel('z (m)'); ylabel('V/m') 

Ppl(l:nkzp-2)=abs(Ezf(nkzp+1:nz,ml)); Ppl(nkzp- 
l:nz)=abs(Ezf(l:nkzp,ml)); 

Pp2 (l:nkzp-2)=eta*abs(Hpf(nkzp+1:nz,ml)); % proper fftshift operation 
Pp2(nkzp-1:nz)=eta*abs(Hpf(1:nkzp,ml)); 

Pp3(l:nkzp-2)=abs(Erf(nkzp+l:nz,ml)); Pp3(nkzp- 
l;nz)=abs(Erf(l:nkzp,ml)); 

subplot(2,1,2), plot(kz,Ppl,'g',kz,Pp2,':r',kz,Pp3,' c ) 
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title('|E_z| (solid); eta*|H_{phi}| (short dash); lE_{rho}| (long 
dash)') 

xlabel('k_z (rad/m)’); ylabel('IFFT1'); figure(1) 
hcpy=input('Enter 1 for Hard Copy: '); 
if ~isempty(hcpy), if hcpy == 1, print; end; end 

clf reset 

subplot(2,1,1), plot(zee,Pr(:,ml),'g',zee,Pz(:,ml), 'r:’) 
if N > 1, 

title(['Power & PSD; Array N=', int2str(N),'; d=',num2str(d), ... 

'm; L=',nim2str(L),'m; phi=',num2str(phi),'deg; f=',... 
ntim2str (f), 'MHz; rho=',num2str (rpl),'m']) 

else, 

title (['Fields and Spectra for L=',nuia2str (L),'la Dipole at f=',... 
nuia2str (f), 'MHz; rho=',num2str (rpl),'m']) 

end 

xlabelCz (m) '); ylabel ('W/m') 

% proper fftshift operation 

Ppl(1:nkzp-2)=Prf(nkzp+1:nz,ml); Ppl(nkzp-1:nz)=Prf(1; nkzp, ml) ; 

Pp2(l;nkzp-2)=Pzf(nkzp+l:nz,ml); Pp2(nkzp-l:nz)=Pzf(l;nkzp,ml) ; 

subplot(2,1,2), plot(kz,Ppl,'g',kz,Pp2,'r:') 

titleCPlots: P_r (solid) P_z (dash)') 

xlabel('k_z (rad/m)'); ylabel ('PSD'); figured) 

hepy=input('Enter 1 for Hard Copy: '); 

if ~isempty(hcpy), if hcpy == 1, print; end; end 

end 

if nplot==8, 

if isempty(rl), error('Enter Rhol First'); end 
r2=input('Enter Rho2 radius to propagate to: '); 

% Selecting closest grid radius 
[rm, m2]=min(abs(rho-r2)); rp2=rho(m2); 

% Initialize Arrays 

TEz=zeros(nz,l); THp=TEz; Ppl=TEz; Pp2=TEz; 

FE2=zeros(nkzp,1); FH2=FE1; 

% Forming numerator of radial transfer functions 

[H, DH]=dhankel(0,krsup*rp2); FE2(nsup)=H; FH2 (nsup)=DH; 

if ~isempty (nsiib), 

[K, DK]=dmodbes(0,krsub*rp2); FE2(nsub)=K; FH2 (nsub)=DK; end 
while 2 

CEz=input('Enter Optimal Estimator Constant for E_z (CEz < 0 to exit): 
'); 

if CEz < 0, break; end 

CHp=input('Enter Optimal Estimator Constant for H_phi & E_rho: '); 

% Computing Riad's Optimal Deconvolution Propagation Transfer Functions 
TEz (l:nkzp)=FE2.*conj (FED ./(FEl.*conj (FED +CEz*EED; 
THp(l:nkzp)=FH2.*conj(FHl)./(FHl.*conj(FHl) + CHp*EHD; 

% Filling in Negative k_z Spectrum Using Symmetry 
TEz(nkzp+l:nz)=flipud(TEz(2:nkzp-l)); 
THp(nkzp+l:nz)=flipuddHp(2:nkzp-D ); 

% Plotting Magnitudes of Transfer Functions 
% proper fftshift operation 

Ppl(l:nkzp-2)=abs(TEz(nkzp+l:nz)); Ppl(nkzp-1:nz)=abs(TEz(l:nkzp)); 

Pp2(l:nkzp-2)=abs(THp(nkzp+l:nz)); Pp2(nkzp-1:nz)=abs(THp(l:nkzp)) ; 
clf reset; plot(kz,Ppl,'g',kz,Pp2,'r:') 
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title ([’Transfer Functions: lTEz(k_z)| (solid);'... 

’ |THp(k_z)l (dash); Rhol=’,num2str(rpl),... 

'm; Rho2=',num2str(rp2),'iii; CEz=',num2str (CEz), . .. 

'; CHp=',num2str(CHp)]) 

xlabel('k_z (rad/m)'); ylabel('Magnitude'); figure(1) 
hepy=input('Enter 1 for Hard Copy: '); 
if ~isempty(hepy), if hepy == 1/ print; end; end 

% Propagating Spectra and Inverse Transforming Fields 
Ezf2=TEz.*Ezf(:,ml); Hpf2=THp.*Hpf(:,ml); Erf2=THp.*Erf(:,ml); 

% Use fft as inverse to ifft for k_z polarity reversal in transforms 
Ez2=fft(Ezf2)/nz; Hp2=fft(Hpf2)/nz; Er2=fft(Erf2)/nz; 

Pr2=-pi*rp2.*real(Ez2.*conj(Hp2)); % rho-directed power/dz 

Pz2=+pi*rp2.*real(Er2.*conj(Hp2)); % z-directed power/dz 

Prf2=-rp2.*real(Ezf2.*conj(Hpf2))/2;% Power Spectral Densities 
Pzf2=+rp2.*real(Erf2.*conj(Hpf2))/2; 

% Computing RMS Errors in Spatial Fields at Rho2 

DEz=Ez(:,m2)-Ez2; Ezerr=100*sqrt((DEz'*DEz)/(Ez(:,m2)'*Ez(:,m2))); 
DHp=Hp(:,m2)-Hp2; Hperr=100*sqrt((DHp'*DHp)/(Hp(:,m2)'*Hp(:,m2))); 
DEr=Er(:,m2)-Er2; Ererr=100*sqrt((DEr’*DEr)/(Er(:,m2)'*Er(:,m2))); 
DPr=Pr(:,m2)-Pr2; Prerr=100*sqrt((DPr'*DPr)/(Pr(:,m2)'*Pr(:,m2))); 
DPz=Pz(:,m2)-Pz2; Pzerr=100*sqrt((DPz'*DPz)/(Pz(:,m2)'*Pz(:,m2))); 

% Plotting Spatial and Spectral Fields at rho=rp2 
eXf ir0S6t 

subplot(2,1,1), plot(zee,abs(Ez(:,m2)),'g',zee,abs(Ez2),'r:') 
title (['E_z at Rho=',num2str(rp2),'m Propagated from Rho=',— 
num2str(rpl),'m; CEz=', num2str(CEz), '; RMS Error=',... 
num2str(Ezerr),’%']) 
xlabel('z (m)'); ylabel('V/m') 

Ppl(l:nkzp-2)=abs(Ezf(nkzp+l:nz,m2)); Ppl(nkzp- 
1:nz)=abs(Ezf(1:nkzp,m2)); 

Pp2 (1:nkzp-2)=abs(Ezf2(nkzp+1:nz)); Pp2(nkzp-1:nz)=abs(Ezf2(1:nkzp)) 

subplot(2,1,2), plot(kz,Ppl,'g’,kz,Pp2,'r:') 

title('Plot Key: Exact (solid) Propagated (dash)') 

xlabel ('k _2 (rad/m)'); ylabel (' I FFT I '); figured) 

hcpy=input('Enter 1 for Hard Copy: d; 

if ~isempty(hepy), if hepy == 1, print; end; end 

eXf 3r6s©t 

subplot(2,1,1), plot(zee,abs(Hp(:,m2)),'g',zee,abs(Hp2),'r:') 
title(['H_{phi} at Rho=',num2str(rp2),'m Propagated from Rho=',... 
num2str(rpl),'m; CHp=',num2str(CHp),'; RMS Error=',... 
num2str(Hperr),'%']) 
xlabel('z (m)'); ylabel('V/m') 

Ppl(l:nkzp-2)=abs(Hpf(nkzp+1:nz,m2)); Ppl(nkzp- 
l:nz)=abs(Hpf(l;nkzp,m2)); 

Pp2(1;nkzp-2)=abs(Hpf2(nkzp+1:nz)); Pp2(nkzp-1:nz)=abs(Hpf2(1:nkzp)) 

subplot(2,1,2), plot(kz,Ppl,'g’,kz,Pp2,'r:') 

title('Plot Key: Exact (solid) Propagated (dash)') 

xlabel('k_z (rad/m)’); ylabel('IFFTl'); figure(1) 

hcpy=input('Enter 1 for Hard Copy: d; 

if ~isempty(hepy), if hepy == 1, print; end; end 

eXf 3r€S6t 

subplot(2,1,1), plot(zee,abs(Er(:,m2)),'g',zee,abs(Er2),'r:') 
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title{['E_{rho} at Rho=',num2str{rp2),’m Propagated from Rho=',... 
num2str(rpl),'m; CEr=CHp=',num2str(CHp ), RMS Error=',... 
num2str(Ererr),'%']) 
xlabel('z (m)'); ylabel(’V/m’) 

Ppl(1:nkzp-2)=abs(Erf(nkzp+1:nz,m2)); Ppl(nkzp- 
1:nz)=abs(Erf(1:nkzp,m2)); 

Pp2(1:nkzp-2)=abs(Erf2(nkzp+1:nz)); Pp2(nkzp-l;nz)=abs(Erf2(l:nkzp)) 

subplot(2,1,2), plot(kz,Ppl,'g',kz,Pp2,'r:') 

title('Plot Key: Exact (solid) Propagated (dash)') 

xlabel('k_z (rad/m)'); ylabel('IFFT|') ; figure(1) 

hcpy=input('Enter 1 for Hard Copy: '); 

if ~isempty(hcpy), if hcpy == 1, print; end; end 

clf reset 

subplot(2,1,1), plot(zee,Pr(:,m2),'g', zee,Pr2,'r:') 

title(['P_{rho} at Rho=',num2str(rp2),'m Propagated from Rho=',... 

num2str(rpl),'m; CEz=',num2str(CEz),'; CHp=',num2str(CHp)__ 

'; RMS Error=',num2str(Prerr),'%']) 
xlabel Cz (m) ') ; ylabel ('W/m') 

Ppl(1:nkzp-2)=Prf(nkzp+1:nz,m2); Ppl(nkzp-l:nz)=Prf(l:nkzp,m2); 

Pp2(1:nkzp-2)=Prf2(nkzp+l:nz); Pp2(nkzp-l:nz)=Prf2(l:nkzp); 

subplot(2,1,2), plot(kz,Ppl,'g',kz,Pp2,'r:') 
titleCPlot Key: Exact (solid) Propagated (dash)') 

xlabel('k_z (rad/m)'); ylabel('P_{rho} PSD'); figure(1) 
hcpy=input('Enter 1 for Hard Copy: '); 
if -isempty(hcpy), if hcpy == 1, print; end; end 

clf reset 

subplot(2,1,1), plot(zee,Pz(:,m2),'g',zee,Pz2,'r:') 
title(['P_z at Rho=',num2str(rp2),'m Propagated from Rho=',... 
num2str ( 2 :pl),'m; CEr=CHp=',num2str(CHp), ... 

'; RMS Error=',num2str(Pzerr),'%']) 
xlabel('z (m) '); ylabel('W/m') 

Ppl(l:nkzp-2)=Pzf(nkzp+l:nz,m2); Ppl(nkzp-1:nz)=Pzf(l;nkzp,m2); 

Pp2(1:nkzp-2)=Pzf2(nkzp+l:nz); Pp2(nkzp-l:nz)=Pzf2(l:nkzp); 

subplot(2,1,2), plot(kz,Ppl,'g',kz,Pp2,'r:') 

titleCPlot Key: Exact (solid) Propagated (dash)') 

xlabel('k_z (rad/m)'); ylabel('P_z PSD'); figure(1) 

hcpy=input('Enter 1 for Hard Copy: '); 

if -isempty(hcpy), if hcpy == 1, print; end; end 

end; end 

if nplot==9, 

if isempty(rl), error('Enter Rhol First'); end 

r2=input('Enter Rho2 for Endpoint of Propagation Error Calculation: '); 

if r2==rl, error('Enter Rho2 Different Then Rhol'); end 

% Selecting closest grid radius 

[rm, m2]=min(abs(rho-r2)); 

if m2 > ml, M=ml:m2; else, M=m2:ml; end 

rp=rho(M); nrp=length(rp); 

FE2=zeros(nkzp,nrp); FH2=FE2; TEz=zeros(nz,nrp); THp=TEz; 

% Forming denominator of radial transfer functions 
for m=l:nrp; 

[H, DH]=dhankel(0, krsup*rp(m)); 

FE2(l:nkzsup,m)=H; FH2(1:nkzsup,m)=DH; 
if -isempty(krsub), % subluminal spectra 

[K, DK]=dmodbes(0,krsub*rp(m)); 
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FE2(nkzsup+l:nkzp,m)=K; 

FH2 (nkzsup+1: nkzp, m) =DK; 
end; end 

FEE=FEl*ones{l,nrp); FHH=FHl*ones(1,nrp); Ppl=zeros(nz,nrp); Pp2=Ppl; 
while 2 

CEz=input('Enter Optimal Estimator Constant for E_z (CEz < 0 to exit): 
'); 

if CEz < 0, break; end 

CHp=input{'Enter Optimal Estimator Constant for H_phi & E_rho: '); 

% Computing Riad's Optimal Deconvolution Propagation Transfer Functions 
TEz(l:nkzp,:)=FE2.*conj(FEE)./(FEE.*conj(FEE) + CEz*EEl); 

THp(l:nkzp,:)=FH2.*conj(FHH)./(FHH.*conj(FHH) + CHp*EHl); 

% Filling in Negative k_z Spectrum Using Symmetry 
TEz(nkzp+1:nz,:)=flipud(TEz(2:nkzp-l,:)); 

THp(nkzp+1:nz,:)=flipud(THp(2:nkzp-l,:)); 

% Plotting Magnitudes of Transfer Functions 

% proper fftshift operation 

Ppl(l:nkzp-2,:)=abs(TEz(nkzp+1:nz,:)); 

Ppl(nkzp-l:nz,:)=abs(TEz(l:nkzp,:)); 

Pp2 U:nkzp-2,:)=abs(THp(nkzp+1:nz,:)); 

Pp2(nkzp-1:nz,:)=abs(THp(1:nkzp, :)); 

clf reset; mesh(rp,kz,Ppl) 

title(t'Deconvolution Function: 1TEz(kz,rho)| with 

'CEz=',num2str(CEz),'; Propagate from Rhol=',nirai2str(rpl),'m']) 
ylabel('k _2 (rad/m)'); xlabel('rho (m)'); view(135,40); figure(1) 
hcpy=input('Enter 1 for Hard Copy: '); 
if ~isempty(hcpy), if hcpy == 1, print; end; end 

clf reset; mesh(rp,kz,Pp2) 

title(['Deconvolution Transfer Function: |THp(kz,rho)| for '... 

'CHp=',num2str(CHp),'; Propagate from Rhol=',num2str(rpl),'m']) 
ylabel(’k_z (rad/m)'); zlabeK'rho (m) ') ; view(135,40) ; figured) 
hcpy=input('Enter 1 for Hard Copy: '); 
if ~isempty(hcpy), if hcpy == 1, print; end; end 

% Propagating Spectra and Inverse Transforming Fields 
Ezf2=TEz.*(Ezf(:,ml)*ones(l,nrp)); 

Hpf2=THp.*(Hpf(:,ml)*ones(l,nrp)); 

Erf2=THp.*(Erf(:,ml)*ones(l,nrp)); 

Ez2=fft(Ezf2)/nz; Hp2=fft{Hpf2)/nz; Er2=fft(Erf2)/nz; 

Pr2=-pi*(ones(nz,1)*rp).*real(Ez2.*conj(Hp2)); % rho-directed power/dz 

Pz2=+pi*(ones(nz,1)*rp).*real(Er2.*conj(Hp2)); % z-directed power/dz 

% Computing RMS Errors in Spatial Fields at Rho2 
DEz=Ez(:,M)-Ez2; DHp=Hp(:,M)-Hp2; DEr=Er(:,M)-Er2; 

DPr=Pr(:,M)-Pr2; DPz=Pz(:,M)-Pz2; 
for m=l:nrp; 
m2=M(m); 

Ezerr (m) =100*sqrt ((DEz (: ,m) '*DEz (: ,m) ) / (Ez (: ,m2) ' *Ez (: ,m2)) ) ; 

Hperr (m) =100*sgrt ((DHp (: ,m) '*DHp (; ,m) ) / (Hpd ,m2) ' *Hp(: ,m2))) ; 

Ererr(m)=100*sqrt((DEr(:,m)'*DEr d/m))/(Er(:,m2)'*Er(:,m2))); 

Prerr(m)=100*sqrt((DPr(:,m)'*DPr(:,m))/(Pr(:,m2)'*Pr(:,m2))); 

Pzerr(m)=100*sqrt((DPz(:,m)'*DPz(: ,m))/(Pz(:,m2)'*Pz(:,m2))); 
end 
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clf reset; plot(rp,Ezerr,'g',rp,Hperr,'r:',rp,Ererr,'c—') 
xlabelCRho (m) '); ylabelCRMS Error (%)'); 
v= axis; v(3)=0; axis (v); figured) 

title(['Error in E_z (solid), H_{phi} (short dash) and E_{rho}'. 
' (long dash) Propagated from Rhol=',num2str(rpl),... 

';CEz=',num2str(CEz),'; CHp=',num2str(CHp)]) 
hcpy=input('Enter 1 for Hard Copy: '); 
if ~isempty(hcpy), if hcpy == 1, print; end; end 

clf reset; plot{rp,Prerr,'g',rp,Pzerr,'r:') 

XlabelCRho (m) '); ylabel ('RMS Error (%)') 
v= axis; v(3)=0; axis(v); figure(1) 

title(['Error in P_r (solid) and P_z (dash) Propagated from Rhol 
num2str(rpl) , ';CEz=',num2str(CEz),'; CHp=',num2str(CHp)]) 
hcpy=input('Enter 1 for Hard Copy: '); 
if ~isempty(hcpy), if hcpy == 1, print; end; end 

end; end; end 
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% EHFLDTST9.M tests dRhoEHFld function routine by comparing computation 
% for sinusoidal I{z)=Im sin k(h-|z|) at selected cylindrical^ 

% radius. (Dipole at the origin with sinusoidal Current Distribution). 

% Using exact integration results for this assumed current 
% given in Jordan & Balmain pp 333-337 and Elliot pp 329-332 for 
% the Exact Fields at radius 1 and radius 2 (> radius 1) and the 
% Computed Ez and Hp fields at radius 2 given the fields at radius 1. 

% Adapted from EHFLDTST.M by M.A. Morgan 6/6/97 
% by D.G. Steenman 12/28/98 

% Mod 9 of EHFLDTST: Calls dRhoEHFld4 which uses different field 
% calculations. 

% Calls dRhoEHFld4 (which calls intphi3 and phisquare) 
clear all 

L=input('Enter dipole length L (meters): '); h=L/2; 

f=input('Enter frequency (MHz): '); 

wl=300/f; k=2*pi/wl; kh=k*h; Im=l; ckh=cos(kh); 

Invk = 1/k; 

disp(['k*h= ',num2str(kh)]); 
disp(['l/k= ',num2str(Invk)]); 

while 1 


zmax=input ('Enter zmax (m) for z—[-zmax ... zmax] plots (0 to stop): ')/ 
if -isempty(zmax), if zmax == 0, break; end; end 
M=input('Enter # z-points for [-zmax ... zmax]: '); 
dz=2*zmax/(M-1); z=(-zmax:dz:zmax); 

rhol=input('Enter rhol (m) for plot (0 to stop): '); 
if ~isempty(rhol)/ if rhol == 0, break; end; end 

rho2min=input {'Enter rho2 min (m) for [rho2min ... rho2max] (0 to stop) : 

’) ; 

if ~isempty (rho2min), if rho2min == 0/ break; end; end 

rho2max=input ('Enter rho2 max (m) for [rho2min ... rho2max] (0 to stop): 

')! 

if ~isempty(rho2max), if rho2max == 0/ break; end; end 

Nrho = input(' Enter # rho2-points for [rho2min ... rho2max]: '); 

dr= (rho2max-rho2min) / (Nrho-1) ; DRHO= (rho2min: dr: rho2max) ; 

NpCalc = ceil(pi*2*rhol/dz) 

Nplow=input('Enter Min number of Phi segments around a thin ring (radius 
1 ) : '); 

if NpCalONplow, 

Np=NpCalc 

else, 

Np=Nplow 

end 

% Exact Field Calculation at Radius 1 
zl=z-h; z2=z+h; 

r=sqrt(rhol*rhol+z.*z); E=exp(-j*k*r); 

Rl=sqrt(rhol*rhol+zl.*zl); El=exp(-j*k*Rl); 

R2=sqrt(rhol*rhol+z2.*z2); E2=exp(-j*k*R2); 
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Ezl=-j*30*Im*(El./Rl+E2./R2-2*ckh*E./r); 

Hpl= (j*Iia/ (4*pi)) * (El+E2-2*ckh*E) /rhol; 

Jz=Hpl; Mp=Ezl; 

% Set Up Loop to do sets of rho2 Field Calculations 
EZRH02 = zeros(length(z),length(DRHO)); ■ 

HPRH02 = zeros(length(z),length(DRHO)); 

EZRH02N = zeros(length(z),length(DRHO)); 

HPRH02N = zeros(length(z),length(DRHO) ) ; 

EZRH02NT = zeros(length(z),length(DRHO)) ; 

HPRH02NT = zeros(length(z),length(DRHO)) ; 

for a = 1:(length(DRHO)), 
rho2 = DRHO(a); 

% Exact Field Calculation at Radius 2 
rrho2=sqrt(rho2*rho2+z.*z); Erho2=exp(-j*k*rrho2) ; 

Rlrho2=sqrt(rho2*rho2+zl.*zl); Elrho2=exp(-j*k*Rlrho2) ; 

R2rho2=sqrt(rho2*rho2+z2.*z2)/ E2rho2=exp(-j*k*R2rho2); 

Ez 2=-j*30*Ini* (Elrho2./Rlrho2+E2rho2. /R2rho2-2*ckh*Erho2. /rrho2); 
Hp2=(j*Im/(4*pi))*(Elrho2+E2rho2-2*ckh*Erho2)/rho2; 

EZRH02(:,a)=Ez2.’; 

HPRH02 (:,a)=Hp2.*; 

% Numerical Integration using dRhoEHFld function for Radius 2 
[Ez2N, Ez2NT, Hp2N, Hp2NT]=dREHF4(Np,k,rhol,rho2,z,z,z,z,l,l,Jz,Mp); 
EZRH02N(:,a)=Ez2N. 

HPRH02N(:,a)=Hp2N.’; 

EZRH02NT(:,a)=Ez2NT.'; 

HPRH02NT(:,a)=Hp2NT.'; 

%clf reset; 
figure(2*a-l) 

plot(z, real(Ez2),'g',z,real (Ez2N),';r',z,real(Ez2NT), ’-.y', ... 
z,imag(Ez2),'g',z,imag(Ez2N),’:r',z,imag(Ez2NT), ... 

z,abs(Ez2), 'g',z,abs(Ez2N),':r',z,abs(Ez2NT),'-.y') 
title(['E_z: Exact (solid); Integ (dot); NullTest (dashdot) at 
{\rho}2=', num2str(rho2),.-. 

' m from {\rho}l= ’ ,num2str (rhol), ' m for L=',num2str (L), ’ m; 
N_z=',num2str(M),... 

'; N_{\phi}=',int2str(Np),f=',num2str(f), ' MHz']); 
xlabel(['z (m), Density N_z = ',num2str(M/(2*zmax)) ]); ylabeK'Real and 
Imag E_z (V/m) '); 
grid; orient landscape; 

%hcpy=input('Print Hard Copy ? (Y/N): ','s'); 

%if ~isempty(hcpy), if hcpy=='Y'Ihcpy=='y',orient landscape; print; end 
end 

%clf reset; 
figure(2*a) 

plot(z,real(Hp2),'g',z,real (Hp2N),’:r',z,real(Hp2NT),'-.y',... 
z,imag(Hp2), ' g', z, imag(Hp2N), ' ;r', z, imag(Hp2NT), '-.y', ... 
z,abs (Hp2), 'g',z,abs (Hp2N), ':r', z,abs (Hp2NT), '-.y') 
title(['H_{phi}: Exact (solid); Integ (dot); NullTest (dashdot) at 
{\rho}2=',num2str(rho2),... 
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' m from {\rho}l= ',num2str (rhol), ' m for L=',niim2str (L), ' m; 

N_ 2 =',num2str(M),... 

N_{\phi}=',int2str(Np),'; f=',num2str(f)/' MHz']); 
xlabel(['z (m), Density N_z = ' ,num2str (M/(2*zmax)) ]) ; ylabeK'Real and 
Imag H_{phi} (A/m)'); 
grid; orient landscape; 

%hcpy=input('Print Hard Copy ? (Y/N); ','s'); 

%if ~isempty(hcpy), if hcpy=='Y'|hcpy=='y', orient landscape; print; 
end; end 

end 

d = 2*length(DRH0); 

%d=0; 

[NRH02G,ZEEG]=iaeshgrid{DRH0./rhol, z); 

RErMagH = ((abs(HPRH02-HPRH02N)))./((abs(HPRH02))) ; 

RErMagE = ((abs(EZRH02-EZRH02N)))./((abs(EZRH02))); 

LSErMagH = sqrt((sum( (abs(HPRH02- 

HPRH02N)) .-'2)) ./ (sum( (abs (HPRH02)) .''2))) *100; 

LSErMagE = sqrt ((Siam Uabs (EZRH02- 

EZRH02N)) .'^2 )) ./ (sum( (abs (EZRH02)) .^2))) *100; 

RErMagHNT = ((abs(HPRH02-HPRH02NT)))./((abs(HPRH02))); 

RErMagENT = {(abs(EZRH02-EZRH02NT)))./((abs(EZRH02))); 

LSErMagHNT = sqrt({sum((abs(HPRH02- 
HPRH02NT)) .''2)) ./ (sum( (abs (HPRH02) ) .''2))) *100; 

LSErMagENT = sqrt( (sum((abs(EZRH02- 
EZRH02NT )) .'■2)) ./ (sum((abs (EZRH02) ) .''2)) ) *100; 


figure(d+1) 

surf(NRH02G,ZEEG,RErMagH) 

title([' Relative Error in Exact to Num Soln for Mag(H_{\phi} ({\rho} 
2)); {\rho} 1= ',num2str(rhol), ' (m) ']) 
xlabeKC {\rho} 2/{\rho} 1 ; N_{\rho}= ' ,num2str (Np) ]) ; 
ylabel([' z (m) ; N_z = ',num2str(M), ' Density N_z = 

' ,num2str (M/ (2*zmax)) ]); 
orient landscape;rotateSd; 

figure(d+2) 

surf(NRH02G,ZEEG, RErMagHNT) 

title([’ Relative Error in Exact to Num Null Test Soln for Mag(H_{\phi} 
({\rho} 2)); {\rho} 1 = ',num2str(rhol), ' (m) ']) 

xlabeKC {\rho} 2/{\rho} 1 ; N_{\rho}= ' ,num2str (Np) ]) ; 
ylabel([' z (m) ; N_z = ',num2str(M),' Density N_z = 

' ,num2str(M/(2*zmax))]); 
orient landscape;rotateSd; 

figure(d+5) 

surf(NRH02G,ZEEG,RErMagE) 

title([' Relative Error in Exact to Num Soln for Mag(E_z ({\rho} 2)); 
{\rho} 1 = ',num2str(rhol),' (m) ']) 

xlabel([' {\rho} 2/{\rho} 1 ; N_{\rho}= ',num2str(Np)]); 
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ylabel([’ z (m) ; N_z = ' ,nirai2str (M), ' Density N_z = 

’,num2str(M/{2*zmax))]); 
orient landscape; rotateSd; 

figure(d+6) 

surf(NRH02G,ZEEG,RErMagENT) 

title([' Relative Error in Exact to Num Null Test Soln for Mag(E_z 
({\rho} 2)); {\rho} 1 = ’ ,nuia2str (rhol), ' (m) ']) 

xlabel{[' {\rho} 2/{\rho} 1 ; N_{\rho}= ' ,nm2str (Np) ]) ; 
ylabel(t' z (m) ; N_z = ' ,nxua2str (M), ' Density N_z = 

' ,nuin2str (M/ (2*zmax)) ]); 
orient landscape; rotateSd; 

figure(d+7) 

plot(DRHO,LSErMagH,'r',DRHO,LSErMagE, ’ — 
g',DRHO,LSErMagHNT,':y',DRHO,LSErMagENT,'-.y') 

title([' Least Squares Error for Mag(H_{\phi}) and Mag(E_z) ; {\rho} 1 
',num2str(rhol),' m ']) 

xlabel([’ {\rho} 2 ; N_{\rho}= ',num2str(Np),' Density N_z = 

' ,nuiti2str (M/ (2*zmax) ) ]); 

ylabel{’ Least Squares Error (%); MagH (solid), MagE (dash), MagH 
NullTest (dots), MagE NullTest (dash-dot) '); 
grid; orient landscape; 


end 


96 




function[EzJ, EzM, HpJ, HpM, Tezj, Tezm, Thpj, Thpm, zOffset] = 
dREHJMA(Np/k,rhol,rho2,zl,z2/Jz,Mp) 

% Method A for Numerically Calculating the Fields at Rho2* 

% Calculate the EzJ, EzM, HpJ, and HpM fields on a cylindrical surface 
of 

% radius 2 due to the equivilent Electric and Magnetic Currents 
% on a cylindrical surface of radius 1. (Radius 2 > Radius 1), 

% Fields are calculated by integrating a constant magnitude current 
% around a thin ring of radius 1 and summing all the contributing 
% pieces of fields from each element in z. 

% The Transforms are an interim step to finding the fields. 

% Includes ability to calculate fields on a z2 smaller than zl. 

% Adapted from M.A. Morgan notes 10 Dec 98 

% Calls intphiS and phisquare2 

% Calculate the Phi Integrations for Nphi, zd, Rhol, and Rho2 
zd = zl-zl(1); 

[u2,u3,u4,u5,v2,v3,v4,v5,dv3,dv4,dv5,w3, w4,w5] = 
intphi3 (Np, k, rhol, rho2, zd) ; 

% Set Up Integration Matrices from Calculated Nphi segmented phi 

% integrations around a thin ring of radius rhol 

znl=length(zl); 

zn2=length(z2); 

izmax = find(zl>(max (z2)))/ 

zOffset=length(zl(izmax)); 

U2= phisquare2(u2,znl,zn2,zOffset); 

U3= phisquare2(u3,znl,zn2,zOffset); 

U4= phisquare2(u4,znl,zn2,zOffset); 

U5= phisquare2(u5,znl,zn2,zOffset); 

V2= phisquare2(v2,znl,zn2,zOffset); 

V3= phisquare2(v3,znl,zn2,zOffset); 

V4= phisquare2(v4,znl,zn2,zOffset); 

V5= phisquare2(v5,znl,zn2,zOffset); 

DV3= phisquare2(dv3,znl,zn2,zOffset); 

DV4= phisquare2(dv4,znl,zn2,zOffset); 

DV5= phisquare2(dv5,znl,zn2,zOffset); 

W3= phisquare2(w3,znl,zn2,zOffset); 

W4= phisquare2(w4,znl,zn2,zOffset); 

W5= phisquare2(w5,znl,zn2,zOffset); 


% Calculate the Components due to the Different Current Components 
% Method A for solving the fields 

dz = abs(zl(2)“Zl (1)); 

% Hp due to Jz (Hp at rho2 due to Hp at rhol) 

hpj = (rho2*U3 - rhol*V3) + j*k*(rho2*U2 ~ rhol*V2); 

Thpj = (rhol*dz/(4*pi))*hpj; 

Hpj = Jz*Thpj; 
clear hpj/ 

% Ez due to Jz (Ez at rho2 due to Hp at rhol) 
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Rezj = {{2+k*k*rho2*rho2)*U3 - 3*rho2*rho2*U5 - 
((rhol/rho2)+2*k*k*rhol*rho2)*V3 ,•. 

+ 6*rhol*rho2*V5 + k*k*rhol*rhol*W3 - 3*rhol*rhol*W5); 

Imezj = j*k*(2*U2 - 3*rho2*rho2*U4 -(rhol/rho2)*V2 + 6*rhol*rho2*V4 - 
3*rhol*rhol*W4); 

Tezj = (rhol*dz*30/(j*k))*(Rezj + Imezj); 

EzJ = Jz*Tezj; 
clear Rezj Imezj; 

% Hp due to Mp (Hp at rho2 due to Ez at rhol) 

Rhpm = (3*rhol*rho2*U5 - k*k*rhol*rho2*U3 
+(2+k*k*rhol*rhol+k*k*rho2*rho2)*V3 ... 

+ k*k*DV3 - 3*(rhol*rhol+rho2*rho2)*V5 -3*DV5 -k*k*rhol*rho2*W3 + 
3*rhol*rho2*W5); 

Imhpm = j*k*(3*rhol*rho2*U4 + 2*V2 ~ 3*(rhol*rhol+rho2*rho2)*V4 - 3*DV4 
+ 3*rhol*rho2*W4); 

Thpm = (rhol*dz/(j*k*4*120*pi*pi))*(Rhpm+Imhpm); 

HpM = Mp*Thpm; 
clear Rhpm Imhpm; 

% Ez due to Mp (Ez at rho2 due to Ez at rhol) 

ezm = (rho2*V3 - rhol*U3) + j*k*(rho2*V2 - rhol*U2); 

Tezm = (rhol*dz/(4*pi))*ezm; 

EzM = Mp*Tezm; 
clear ezm; 
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function [U1,U2,U3,U4,U5,DU3,DU4,DU5,V1,V2,V3,W3,W4,W5] = 
intphi4 (N, k, rhol, rho2, z2) 

% Program uses trapezoidal rule to evaluate integrals on phi=0 to 2pi 
circular contour 

% of radius rhol in the x-y plane with field point at (rho2, z2,phi=0) of 
even functions: 

% un(phi)= (1/R^n) exp(-jkR) 

% vn(phi)= cos(phi) un(phi) 

% Dun(phi)= (z2''2) un(phi) 

% and wn(phi)= cos"'2{phi) un(phi) 

% where R=sqrt {rhol''2+rho2''2-2rhol*rho2*cos (phi)+z2''2) 

% Adapted from IntPhi.m and IntPhi2.m by M.A. Morgan 12/28/98 

dp=pi/(N-1); phi=0:dp:pi; % double the integral from 0 to pi 
a=rhol*rhol+rho2*rho2+z2.''2; b=2*rhol*rho2; cp=cos(phi); 

A=a. '*ones (1, length (cp)) ; D= (z2 .''2) . '*ones (1, length (cp)); 

CP=ones(length(z2),1)*cp; 

R=sqrt(A-b*CP) ; 
ul=exp(-j*k*R)./(R); 

u2= ul./R; u3= u2./R; u4= u3./R; u5= u4./R; 
vl=CP.*ul; v2=CP.*u2; v3=CP.*u3; v4=CP.*u4; v5=CP.*u5; 
w3=CP.*v3; w4=CP.*v4; w5=CP.*v5; 
du3=D.*u3; du4=D.*u4; du5=D.*u5; 
c=2*dp; ds=c*ones(N,1); ds(l)=c/2; ds(N)=c/2; 

Ul=ul*ds; U2=u2*ds; U3=u3*ds; U4=u4*ds; U5=u5*ds; 

Vl=vl*ds; V2=v2*ds; V3=v3*ds; 

W3=w3*ds; W4=w4*ds; W5=w5*ds; 

DU3=du3*ds; DU4=du4*ds; DU5=du5*ds; 
clear A D CP R u* v* w* dv* ds cp phi; 
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function A = phisquare4 (Vecn,Nzl,N22,MC,N, zOffsetMax, zOffsetMin) 

% 

% Mod4: (length of Intdzl) *N== (length of d22) : Also, .use MATLAB 
% ’toeplitz’^ function instead of spdiags, 

% 

% Mod3: Takes different dzl | dz2 sizes AND 
% different zlMax | z2Max values 

% into account 

% 

% PhisquareS takes a vector un and makes a rectangular diagonal matrix 
% out of its elements. 'Vecnd) ’ is the main diagonal, 'Vecn(2) ' 

% is the first diagonal off the main, both above and below 
% the diagonal. 'Veen* is a column vector. After making a 
% large rectangular matrix, specific rows are choosen to account 
% for the difference in dzl compared to dz2. 

% If they are the dz's are equal, this program makes a banded 
% Toeplitz matrix. 

% by D.G. Steenman 12/30/98 


% Set up elements for each row 

U1 = Veen.»; 

ANz = length(Veen); 

U = toeplitz(U1,U1); 

% First, Deal with Offset 

% Second, Choose Selected Columns that correspond to 2R2 points 
% Third, reshape and add N*MC rows to average the integrations 
% Fourth, squeeze and transpose to correct orientation for A 

if (zOffsetMax<=0), 

U=U(:, (~zOffsetMax+1) : (ANz+zOffsetMin)); 
else, 

U=U( (zOffsetMax+1) ; (ANz-zOffsetMin), ;) ; 

end 


U=U(:,((0:Nz2-l)*N + 1)); 
U=U.^; 

U=reshape(U,Nz2,N*MC,Nzl); 
U=sum(U, 2) ; 

Unsqueeze (U); 

A=U.»; 

clear U* ANz; 
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